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REPORT  SUMMARY  AND  DISCUSSION 


Electromagnetic  penetration  and  scattering  problems  are  difficult 
to  treat  with  many  analytical  or  numerical  methods  because  of  the  Ina¬ 
bility  of  these  methods  to  simply  deal  with  the  effects  of  structure 
materials,  apertures,  curvatures,  corners,  and  Internal  contents.  In  two 
earlier  RADC  contracts,  F30602-77-C-0163  and  F30602-79-C-0039,  IIT  Research 
Institute  (IITRI)  Investigated  the  application  of  a  new  approach  for  the 
direct  modeling  of  very  complex  electromagnetic  Interaction  problems:  the 
finite-difference,  time-domain  (FD-TD)  solution  of  Maxwell's  equations. 

The  FD-TD  method  has  key  advantages  relative  to  available  modeling  approaches. 
These  advantages  permit  It  to  accurately  treat  complex  problems  that  are 
beyond  the  scope  of  solution  by  any  other  method.  The  ultimate  aim  of 
research  In  this  area  Is  to  develop  an  accurate,  easily-used,  general  com¬ 
puter  program  solving  for  either  electromagnetic  field  penetration,  scat¬ 
tering,  or  radiation  for  arbitrary  metal /dielectric  structures  spanning  up 
to  10  or  more  wavelengths  In  three  dimensions  with  a  spatial  resolution 
better  than  0.1  wavelength. 

In  order  to  more  fully  determine  the  usefulness  of  the  FD-TD  method, 

RADC  thought  it  Is  desirable  to  distribute  this  technique  to  as  wide  a 
range  of  users  as  possible  so  that  It  can  be  tested  by  actual  Implementation. 
The  overall  objectives  of  algorithm  development  In  this  case  are  to  allow 
RADC  to  write  a  user-oriented  computer  program  for  the  FD-TD  technique. 

The  goals  of  the  present  IITRI  research  effort  for  RADC,  Contract 
F30602-80-C-0302,  Included  the  development  of  specific  algorithms  of  high 
Importance  to  help  provide  a  flexible,  s1mple-to-use,  and  highly  accurate 
user-oriented  FD-TD  computer  program.  To  meet  these  goals,  IITRI  tested 
five  key  Improvements  In  the  FD-TD  algorithm  during  this  effort,  including 
the  following: 


1 .  Total “ field/ sea tterect-fi eld  lattice  division. 

This  permits  a  very  high  computational  dynamic  range  to  accurately 
model  fields  within  shadow  zones  or  cavities.  This  further  permits 
programming  of  variable  angle  of  incidence  and  the  second-order 
correct  radiation  condition,  summarized  helow. 

2 .  Variable  angle  of  incidence. 

For  two-dimensional  problems,  this  permits  a  single  data  card  or 
Fortran  statement  to  specify  in  a  very  accurate  manner  the  angle 
of  incidence  of  a  plane  wave  illuminating  a  structure.  For  three- 
dimensional  problems,  both  the  angle  of  incidence  and  polarization 
could  be  specified.  There  is  np_  requirement  to  rotate  the  geom¬ 
etry  of  the  interacting  structure  in  the  FD-TD  lattice. 

3.  Second-order  accurate  radiation  condition. 

This  reduces  the  uncertainty  of  the  final  computed  results  by  as 
much  as  ten-to-one.  FD-TD  computations  using  this  radiation  con¬ 
dition  now  have  estimated  field-magnitude  uncertainties  of  better 
than  12.52  (+0.2  dB)  versus  previous  uncertainties  of  +102-+15% 

(+1  dB). 

4.  Magnitude  and  phase  computation  condition  for  the  sinusoidal 

steady  state _ _ 

This  permits  accurate  determination  of  the  magnitude  and  phase 
of  FD-TD  computed  fields  at  any  desired  points  for  later  use  in 
computations  involving  scatteiing,  radiation,  or  coupling  to  wires. 
This  approach  avoids  any  ambiguity  due  to  either  a  possible  DC 
offset  of  the  fields  or  the  repetitive  nature  of  the  sinusoidal 
waveform. 

5 .  Near-to-far  field  transformation. 

This  permits  the  far  scattered  fields  and  radar  cross  section  of 
arbitrary  structures  modeled  by  the  FD-TD  method  to  be  easily  and 
accurately  determined.  Observed  accuracy  of  the  radar  cross  section 
using  this  feature  is  in  the  order  of  +.12  (,t0.09  dB). 


1v 


These  FD-TD  algorithm  improvements  are  documented  in  this  report.  In  ad¬ 
dition  to  the  algorithm  improvements  tested  by  IITRI,  this  report  also  summa¬ 
rizes  a  FD-TD  feature  which  has  recently  appeared  in  the  literature  that  per¬ 
mits  computation  of  the  coupl  ; "g  oi  ci*.;-ents  to  thin  wires  and  struts. 

The  conclusions  of  this  report  are  as  -follows: 

1.  The  accuracy  of  the  pure  FO-TD  method  for  electromagnetic  interaction 
problems  can  reach  the  high  levels  previously  attained  only  by  method- 
of-moments  (MOM)  approaches  when  the  second-order  accurate  radiation 
condition  is  used  in  the  FD-TD  algorithm.  The  FD-TD  method  retains  its 
significant  advantages  over  MOM  in  terms  of  the  much  larger  electrical 
size  and  greater  complexity  of  the  structures  that  can  be  modeled. 

?.  The  specification  of  variable  angle  of  incidence  and  polarization  of 
an  illuminating  wave  can  be  achieved  with  the  FD-TD  method  using  only 
a  single  data  card  or  Fortran  statement.^ 

The  total-field/scattered-field  regional  division  of  the  FD-TD  lattice 
can  be  successfully  implemented  and  offers'  the  significant  advantage 
of  a  high  computational  dynamic  range.  In  addition,  this  lattice  divi¬ 
sion  provides  a  framework  for  programming  variable  wave  incidence  and 
polarization,  Improved  radiation  conditions,  and  the  near-field  to 
far-field  transformation  for  scattering  problems. 

4.  The  near-to-far  field  transformation  along  a  rectangular  virtual  sur¬ 
face  surrounding  a  scatterer  makes  it  possible  to  use  the  FD-TD  method 
to  compute  the  far  mattered  fields  and  radar' cross  section  of  complex, 
arbitrary  structures  with  great  precision. 

It  is  the  opinion  of  the  authors  of  this  report  that  the  FD-TD  method  deserves 
additional  investigation  to  probe  just  what  are  the  limits  of  application  of 
this  extremely  promising  approach  to  accurately  model  electromagnetic  penetra¬ 
tion,  scattering,  and  radiation  problems. 

This  incident  wave  specification  is  now  as  simple  for  the  FD-TD  metnod  as  it 
has  been  with  MOM.  However,  the  FD-TD  approach  requires  re-running  the  entire 
problem  for  each  new  incident  wave  angle.  With  MOM,  only  a  single  Inversion 
of  the  system  matrix  is  required.  Subsequently,  arbitrary  wave  excitation  is 
treated  as  a  simple  matrix  multiplication  of  the  equivalent  excitation  vec¬ 
tor.  MOM  therefore  permits  a  conceptually  simpler  treatment  of  the  variable 
wave  incidence  problem. 
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PREFACE 


I IT  Research  Institute  (IITRI)  Is  pleased  to  submit  this  Final  Report 
on  "User's  Code  for  the  Finite-Difference  Time-Domain  Method"  to  Rome  Air 
Development  Center  (RADC/RBCT).  The  report  covers  work  performed  by  IITRI 
under  Air  Force  Contract  No,  F30602-80-C-0302,  designated  as  IITRI  Project 
No.  E6502.  This  report  covers  details  of  the  technical  work,  including 
relevant  theory  and  numerical  results.  Appendix  A  provides  a  listing  of 
the  Fortran  computer  program  used  to  obtain  the  results  of  this  report. 

The  principal  investigator  on  this  program  was  Dr.  Allen  Taflove. 

The  co-investigator  was  Dr.  Korada  Umashankar.  The  project  duration  was 
1  September  1980  to  31  August  1981. 

Respectfully  submitted, 

I IT  RESEARCH  INSTITUTE 

Allen  Taflove,  Ph.D. 

Senior  Engineer 
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Air  Force  Contract  Ho.  F30602“80-C-0302 
IITRI  Project  No.  E6502 
1  September  1930  -  31  August  1931 
FINAL  REPORT 


1.0  INTRO  .DUaiON 

Electromagnetic  penetration  and  scattering  problems  are  difficult  to 
treat  with  many  analytical  or  numerical  methods  because  of  the  inability 
of  these  methods  to  simply  deal  with  the  effects  of  structure  materials, 
apertures,  curvatures,  corners,  and  internal  contents.  Usually,  only 
relatively  simple  geometries  are  studied  in  an  attempt  to  gain  Insight  into 
the  key  interaction  mechanisms  and  to  allow  an  indirect  estimate  of  the 
interaction  for  more  complicated  problems. 

In  earlier  RADC  Contracts  F30602-77-C-0163  and  F30602-79-C-0039, 

IITRI  investigated  the  application  of  a  new  approach  for  the  direct  model¬ 
ing  of  very  complex  electromagnetic  interaction  problems:  the  finite- 
difference,  time-domain  (FD-TD)  solution  of  Maxwell's  equations.  The  FD- 
TD  method  treats  the  illumination  of  a  structure  as  an  ini tial -value 
problem.  At  t  =  0,  a  plane-wave  source  of  frequency,  f,  is  assumed  to 
be  turned  on.  The  propagation  of  waves  from  this  source  is  simulated  by 
solving  a  finite-difference  analog  of  the  time-dependent  Maxwell's  equa¬ 
tions  on  a  lattice  of  cells,  Including  the  structure.  Time-stepping  is 
continued  until  the  sinusoidal  steady  state  is  achieved  at  each  cell.  The 
field  envelope,  or  maximum  absolute  value,  during  the  final  wave-cycle  of 
time-stepping  is  observed  to  obtain  the  magnitude  and  phase  of  the  steady- 
state  field  at  any  point. 

This  method  has  two  key  advantages  relative  to  available  modeling 
approaches.  First,  it  is  simple  to  implement  for  complicated  metal/die- 
lectric  structures  because  arbitrary  electrical  parameters  can  be  assigned 
to  each  lattice  cell  using  a  data  card  deck.  Second,  its  computer  memory 
and  running  time  requirement  is  not  prohibitive  for  many  complex  structures 
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of  interest .  In  the  RADC  work,  IITRI  has  shown  the  FD-TD  method  to  be 
capable  of  accuraLely  sclvituj  for  more  than  one  million  unktiowii  field 

components  wlthl^n  a  few  minutes  on  an  array-processing  computer.  Con¬ 
sistently,  a  +  1-dB  accuracy  relative  to  known  analytical  and  experimental 
bench  marks  has  been  achieved  for  a  variety  of  dielectric  and  metal  geom¬ 
etries. 

The  ob.iectivo  of  IITRI's  previous  RADC  studies  was  to  evaluate  the 
suitability  of  the  FD-TD  method  to  determine  the  amount  of  electromagnetic 
coupling  through  an  aperture  into  an  enclosed  conducting  container  and  the 
interaction  and  coupling  of  the  penetrating  fields  with  Internal  electronics. 
Two  specific  container  models  were  used  for  the  evaluation.  The  first,  a 
simple  conducting  cylinder  with  one  open  end.  The  other,  the  complex  guid¬ 
ance  section  of  a  missile.  Each  of  these  two  configurations  was  modeled 
to  calculate  the  electromagnetic  field  coupled  into  the  structure. 

The  ultimate  aim  of  research  in  this  area  is  two-fold.  First,  develop 
an  accurate, easily-used,  general  code  solving  for  either  electromagnetic 
field  penetration,  scattering,  or  radiation  for  arbitrary  metal/dielectrlc 
structures  spanning  up  to  10  or  more  wavelengths  1h  three  dimensions  with 
a  spatial  resolution  better  than  0.1  wavelength.  Second,  develop  a  more 
sophisticated  intuitive  understanding  of  basic  wave  interaction  mechanisms 
in  time  domain,  such  as  transient  propagation  through  beyond-cutoff  cavity 
interiors,  field  build-up  at  edges,  convergence  to  the  sinusoidal  steady 
statft  scattering,  and  radiation. 

In  order  to  more  fully  determine  the  usefulness  of  the  FD-TD  method, 

RADC  thought  that  it  is  desirable  to  distribute  this  technique  to  as  wide 
a  range  of  users  as  possible  so  that  it  can  be  tested  by  actual  implemen¬ 
tation.  However,  the  FD-TD  computer  programs  previously  developed  by  IITRI 
for  RADC  were  research  oriented  rather  than  user  oriented,  i.e.,  they  were 
written  to  apply  to  fairly  specific  types  of  interaction  structures  rather 
than  completely  general  structures.  Further,  these  FD-TD  programs  were 
optimized  for  a  specific  vector  array  processing  computer,  the  Control  Data 


Cyber  203,  in  order  to  miniinize  the  cost  of  purchased  computer  time.  The 
Fortran  used  for  these  programs  was  specialized  to  benefit  from  the  machine- 
specific  features  of  the  Cyber  203,  and  is  not  directly  usable  by  common 
scalar-processing  computers  such  as  the  Honeywell  6000  series,  CDC  6600, 

CDC  Cyber  76,  or  IBM  370  series. 

The  overall  objectives  of  algorithm  development  in  this  case  are  to 
allow  RADC  to  write  a  user-oriented  program  for  the  FD-TD  technique.  Such 
a  program  would  ideally  buffer  the  user  from  the  complexities  Involved  in 
specifying  an  interaction  geometry  in  the  form  needed  by  the  main  program. 
Further,  such  a  program  would  be  suitable  for  arbitrary  shaped  bodies, 
material  parameters,  incident-wave  angle  of  incidence  and  direction  of 
polarization,  and  time  dependence  on  the  incident  wave.  Finally,  such  a 
computer  program  would  be  machine-independent,  i.e.,  written  in  a  standard 
language  such  as  Fortran  IV,  so  that  implementation  on  a  very  wide  variety 
of  computers  would  be  easily  possible.  However,  the  program  would  still 
be  structured  to  make  the  best  economy  of  computer  resources  such  as 
memory  storage  and  program  execution  time. 

The  goals  of  the  present  IITRI  research  effort  for  RADC  include  the 
development  of  specific  algorithms  of  high  importance  to  help  provide  a 
flexible,  simple-to-use,  and  highly  accurate  user-oriented  FD-TD  computer 
program.  IITRI  has  tested  five  key  improvements  in  the  FD-TD  algorithm 
during  this  effort,  and  reports  on  these  developments  in  this  document. 

The  following  is  a  listing  of  these  algorithm  developments,  including 
comments  indicating  the  usefulness  of  each  development. 

1 •  Total-field/scattered-field  lattice  division. 

This  permits  a  very  high  computational  dynamic  range  to  accurately 
model  fields  within  shadow  zones  or  cavities.  This  further  permits 
programming  of  variable  angle  of  incidence  and  the  second-order 
correct  radiation  condition,  summarized  below. 
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2.  Variable  angle  of  Incidence. 

For  two-dimensional  problems,  this  permits  a  single  data  card 
or  Fortran  statement  to  specify  in  a  very  accurate  manner  the 
angle  of  Incidence  of  a  plane  wave  Illuminating  a  structure. 

For  three-dimensional  problems,  both  the  angle  of  incidence  and 
polaritatlon  could  be  specified.  There  is  no  requirement  to 
rotate  the  geonietry  of  the  interacting  structure  in  the  FD-TD 
lattice. 

3 .  Second-ordur  accurate  radiation  condition. 

This  reduces  the  uncertainty  of  the  final  computed  results  by 
as  much  as  ten-to-one.  FD-TD  computations  using  this  radiation 
condition  now  have  estimated  field-magnitude  uncertainties  of 
better  than  +2.5%  (+^0.2  dB)  versus  previous  uncertainties  of 
+10%-+! 5%  (+1  dB). 

4.  Magnitude  and  phase  computation  condition  for  the  sinusoidal 

steady  state _ _ _ 

This  permits  determination  of  the  magnitude  and  phase  of  FD-TD 
computed  fields  at  any  desired  points  for  later  use  In  computations 
Involving  scattering,  radiation,  or  coupling  to  wires.  This 
approach  avoids  any  ambiguity  due  to  either  a  possible  DC  offset 
of  the  fields  or  the  repetitive  nature  of  the  sinusoidal  waveform. 

Near-to-far  field  transformation. 

This  permits  the  far  scattered  fields  and  radar  cross  section  of 
arbitrary  structures  modeled  by  the  FD-TD  method  to  be  easily  and 
accurately  determined.  Observed  accuracy  of  the  radar  cross  section 
using  this  feature  1s  1n  the  order  of  (+0.09  dB). 
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In  addition  to  the  five  algorithm  Improvements  tested  by  IITRI ,  this 
report  will  also  summarize  a  FD-TD  feature  which  has  recently  appeared  in 
the  literature  that  permits  computation  of  the  coupling  of  currents  to 
thin  wires  and  struts. 

Section  2.0  of  this  report  will  provide  a  complete  summary  of  the 
theory  and  basic  algorithms  of  the  FD-TD  method,  including  the  five  new 
features  tested  by  IITRI  and  the  thin-wire  coupling  model.  Sub-sections 
which  contain  discussion  of  these  new  features  will  be  denoted  by  a  double 
asterisk  (**)  appearing  before  the  headings.  Section  3.0  of  this  report 
will  provide  examples  of  computed  results  which  illustrate  each  of  the  key 
features  of  the  FD-TD  method  examined  to  date.  Again,  examples  of  new 
features  will  be  denoted  by  a  double  asterisk  appearing  before  the  appro¬ 
priate  heading.  Appendix  A  provides  a  standard  Fortran  listing  of  a  two- 
dimensional  FD-TD  computer  program  which  illustrates  the  key  features 
tested  by  IITRI  during  this  research  effort. 

2 • 0  THEORY  AND  BASIC  ALGORITHMS  OF  THE  FD-TD  METHOD 

2 • 1  Ideas  Behind  the  FD-TD  Method 

The  FD-TD  method  Is  a  direct  solution  of  Maxwell's  time-dependent 
curl  equations.  As  shown  In  Figure  1,  the  goal  Is  to  model  the  propagation 
of  an  electromagnetic  wave  Into  a  volume  of  space  containing  a  dielectric 
or  conducting  structure.  By  time-stepping,  I.e.,  repeatedly  implementing 
a  finite-difference  analog  of  the  curl  equations  at  each  cell  of  the 
corresponding  space  lattice,  the  Incident  wave  Is  tracked  as  1t  first 
propagates  to  the  structure  and  then  Interacts  with  it  via  surface-current 
excitation,  diffusion,  penetration,  and  diffraction.  Wave-tracking  Is 
completed  when  the  desired  late-tliiie  or  sinusoidal  steady-state  behavior 
is  observed  at  each  lattice  cell.  The  rationale  for  this  procedure  Is 
that  It  achieves  simplification  by  analyzing  the  interaction  of  the  wave- 
front  with  portions  of  the  structure  surface  at  a  given  Instant  In  time, 
rather  than  attempting  a  simultaneous  solution  of  the  entire  problem. 
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Lattice  Truncation  Plone 


Figure  1.  TIME-DOMAIN  WAVE-TRACKING  CONCEPT  OF  THE  FD-TD  METHOD 


T1me*stepp1ng  for  the  FD-TD  method  1s  accomplished  by  an  explicit 
f1nite-d1ffei"ence  procedure  due  to  Yee  [1].  For  a  cubic-cell  space 
lattice,  this  procedure  Involves  positioning  the  components  of  E  and  H 
about  a  unit  cell  of  the  lattice,  as  shown  In  Figure  2,  and  evaluating  E 
and  H  at  alternate  half-time  steps.  In  this  manner,  centered  difference 
expressions  can  be  used  for  both  the  space  and  time  derivatives  to  attain 
second-order  accuracy  In  the  space  and  time  Increments  without  requiring 
simultaneous  equations  to  compute  the  fields  at  the  latest  time  step. 

The  finite-difference  formulation  of  the  FD-TD  method  allows  the 
straightforward  modeling  of  the  surfaces  and  Interiors  of  arbitrary  die¬ 
lectric  or  conducting  structures.  The  structure  of  interest  Is  mapped 
Into  the  space  lattice  by  first  choosing  the  space  Increment  and  then 
employing  a  data  card  deck  to  assign  values  of  permittivity  and  conductivity 
to  each  component  of  E.  No  special  handling  of  electromagnetic  boundary 
conditions  at  media  Interfaces  Is  required  because  the  curl  equations 
generate  these  conditions  in  a  natural  way  by  themselves,  Therefore,  the 
basic  computer  program  need  not  be  modified  to  change  from  structure  to 
structure.  In  this  manner.  Inhomogeneities  or  fine  details  of  the  structure 
can  be  modeled  with  a  maximum  resolution  of  one  unit  cell;  thin  surfaces 
can  be  modeled  as  stepped-edge  sheets.  Figure  3  shows  an  arbitrary  three- 
dimensional  scatterer  embedded  In  a  FD-TD  lattice. 

The  explicit  formulation  of  the  FD-TD  method  Is  particularly  suited  for 
programming  with  minimum  storage  and  execution  time  using  recently  developed 
array-processing  computers.  First,  the  required  computer  storage  and  running 
time  Increases  only  linearly  with  N,  the  total  number  of  unknown  field 
components.  Computer  techniques  (such  as  the  method-of-iiionients)  which  re¬ 
quire  the  solution  of  simultaneous  equations  usually  have  a  storage  require¬ 
ment  proportional  to  and  a  running  time  proportional  to  N‘^-N^  [2]. 

Second,  since  all  FD-TD  operations  are  explicit  and  can  be  performed  in 
parallel,  rapid  array-processing  techniques  can  be  readily  applied.  As 
will  be  demonstrated  later,  these  can  be  employed  to  solve  for  more  than 
10®  field  components  In  a  single  FD-TD  problem,  as  opposed  to  a  maximum  of 
about  10  field  components  for  conventional  approaches  using  simultaneous- 
equation  solutions  [2]. 
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Yl’c  ii|'pli('d  I  h('  I'D- 111  method  Lu  c:oiii|iu  t,c  Uio  wave  form'.  oi  fM  and  Tf 

pulses  scattered  from  infinitely  long,  rectangular  cross  section,  con¬ 
ducting  cylinders  [1],  Other  workers  investigated  electromagnetic-pulse 
(EMP)  interactions  in  time-varying  inhomogeneous  media  [3];  with  metallic 
bodies  of  revolution  ['ll;  with  detailed  models  of  aircraft  [5],  [6]*,  with 
lossy  dielectric  stuctures  [7l;  and  with  thin  struts  and  wires  [8].  Four 
distinct  problems  emerged  in  the  process  of  adapting  the  FD-TO  method  to 
model  realistic  situations: 

1)  Lattice  Truncation  Conditions.  The  field  components  at  the  lattice 
truncation  planes  cannot  be  determined  directly  from  the  Maxwell 's-equations 
analog  and  must  be  computed  using  an  auxiliary  truncation  condition.  However, 
care  must  be  exercised  because  this  condition  must  not  cause  excessive 
spurious  reflection  of  waves  scattered  outward  by  the  structure  modeled. 

The  goal  is  to  formulate  truncation  planes  as  close  as  possible  to  the  struc¬ 
ture  (to  minimize  computer  storage),  and  yet  achieve  virtual  Invisibility 
of  these  planes  to  all  possible  waves  w*''hin  the  lattice. 

2)  Incident  Plane-Wave  Source  Condition.  The  simulation  of  either  a 
plane-wave  pulse  or  single-frequency  plane  wave  at  arbitrary  angles  of 
incidence  should  not  take  excessive  storage  nor  cause  spurious  wave  reflec¬ 
tions.  The  former  would  occur  if  the  incident  wave  is  programmed  as  an 
initial  condition;  the  latter  would  occur  if  the  Incident  wave  is  programmed 
as  a  fixed  field  excitation  along  a  single  lattice  plane. 

3)  Sinusoidal  Steady-State  Information.  Such  data  can  be  obtained 
either  by  a)  directly  programming  a  single-frequency  incident  plane  wave 
or  b)  performing  a  separate  Fourier  transformation  step  on  the  pulse  wave¬ 
form  response.  Both  methods  require  time-stepping  to  a  maximum  time  equal 
to  several  wave  periods  at  the  desired  frequency.  The  second  method  has 
two  additional  requirements.  First,  a  short-rise-time  pulse  suffers  from 
accumulating  waveform  error  due  to  overshoot  and  ringing  as  it  propagates 
through  the  space  lattice.  This  leads  to  a  numerical  noise  component  which 
should  be  filtered  before  Fourier  transformation.  Second,  Fourier  trans- 
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formation  of  many  lattice-cell  field- versus- time  waveforms  (each  probably 
extending  over  many  hundreds  of  time  steps)  would  significantly  add  to 
the  total  requirements  for  computer  storage  and  execution  time. 

4)  Total-Field  versus  Scattered-Field  Formulation.  A  choice  exists 
in  whether  to  finite-difference  only  the  scattered  field  instead  of  the 
total  field  (at  each  lattice  cell).  The  scattered-field  approach  may  lead 
to  a  relatively  superior  lattice  truncation  condition  [4].  However,  the 
total-field  approach  is  more  useful  in  determining  the  fields  penetrating 
structures  having  shielding  properties,  or  the  fields  in  the  shadow  regions 
of  scatterers,  where  the  total  field  can  diminish  to  levels  far  below  the 
incident.  Scattered-field  codes  have  traditionally  run  into  numerical 
"noise"  problems  for  such  cases  since  they  achieve  interior  or  shadow-zone 
field  reduction  by  the  subtraction  of  nearly  equal  scattered  and  incident 
field  quantities.  Computed  shielding  or  shadowing  of  more  than  30  dB  may 
be  difficult  to  achieve  in  this  manner  because  of  a  residual  "noise"  floor 
inherent  in  this  subtraction  process.  A  total-field  approach  does  not 
suffer  from  the  subtraction-noise  problem  and  hence  is  suitable  for  computing 
field  penetration  within  shielded  structures  or  into  shadow  zones. 

Previous  work  by  IITRI  described  initial  efforts  to  solve  the  first 
three  problems  above  for  the  case  of  a  total-field  FD-TD  program  employing 
a  cubic-unit-cell  space  lattice  [9]  -  [13],  Simple  truncation  conditions 
were  developed  for  two  and  three-dimensional  lattices  that  reduced  the  re¬ 
flection  coefficient  of  closely  positioned  truncation  planes  to  the  order 
of  0.1  for  waves  of  arbitrary  Incidence.  A  plane-wave  source  condition  was 
described  that  allowed  generation  of  an  arbitrary  pulsed  or  sinusoidal  wave 
(at  a  fixed  angle  of  incidence)  without  requiring  any  additional  storage 
and  without  causing  spurious  wave  reflections.  Finally,  it  was  shown  that 
sinusoidal  steady-state  data  could  be  computed  using  the  FD-TD  method  by 
directly  programming  a  single-frequency  Incident  plane  wave  and  time-step- 
ping  to  the  steady-state  over  several  periods  of  the  incident  wave.  Ob¬ 
served  accuracy  was  +1  at  points  of  electric-field  maxima,  and  +1  lattice 
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cell  ill  locating  el  oc  tni  c-fi  al  d  maxima  and  iiiinima,  when  tiui  lal.t.  i  1  1 

size  was  chosen  to  be  less  than  0.05  X,  as  discussed  in  detail  in  [12]  and 
[13]. 

The  following  sub-sections  will  provide  a  full  review  of  the  latest 
formulation  of  the  FD-TD  method,  including  the  features  developed  during 
the  current  research  effort  as  well  as  the  thin-wire  coupling  model.  These 
new  developments  will  be  indicated  by  a  double-asterisk  (**)  appearing  before 
the  sub-section  heading. 

2*2  Computational  Details  for  a  Uniform,  Cubic  Lattice 

2.2.1  Basi c_ Jy s_^tejn  of  Equations 

Using  the  MKS  system  of  units,  and  assuming  that  the  material  parameters, 
p,  L  and  a  are  Independent  of  time,  the  following  system  of  scalar  equations 
is  equivalent  to  Maxwell's  equations  in  the  rectangular  coordinate  system 
(x,  y,  z): 
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Yee  [1]  originally  introduced  a  set  of  finite-difference  equations  for 
Lhe  system  of  Lcjuntions  (la)  -  (If).  Following  Yee's  notation,  we  denote  a 
space  point  in  a  cubic  lattice  as 

(i,J,k)  ■•=  (if;,j6,l<iS)  (2) 

and  any  function  of  space  and  time  as 

F'\l,j,k)  -  F(i>i,j6,k6,n6t),  (3) 

where  =  “-x  =  l^y  =  vSz  is  the  space  increment,  (St  is  the  time  Increment,  and 
1,  J,  k,  and  n  are  integers.  Yee  used  centered  finite-difference  expressions 
for  the  space  and  time  derivatives  that  are  both  simply  programmed  and  second- 
order  accurate  in  (S  and  in  (St,  respectively; 


;ix  (S 


+  0(6^) 


J  ,K)_  ,  t  0(J,2,  ,5, 

To  achieve  tlie  accuracy  of  Equation  (4),  and  to  realize  all  of  the  space 
derivatives  of  Equations (la)  -  (If),  Yee  positioned  the  components  of  E  and  H 
about  a  unit  cell  of  the  lattice  as  shown  in  Figure  2.  To  achieve  the  at(  u- 
racy  of  Equation  (5),  he  evaluated  f  and  H  at  alternate  half  time  steps. 

Many  electromagnetic  interaction  problems  involve  nonpermeable  media 
and  can  be  approached  using  a  fixed  time  step  and  space  increment.  For  such 
problems,  the  quantity  dt/ p(  1  ,J ,k)i.'i  is  constant  for  all  (i,J,k)  of  the 
lattice,  and  the  Yoe  system  of  equations  can  be  simplified  to  reduce  computer 
time  in  the  following  manner.  We  define  the  constants: 


R  -  i.St/2c^ 

(6a) 

(6b) 
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%  *  fiVWpS  (6c) 

1  -  Rc(ni)/i.^(m) 

'  r+  Ra(m)/e'^M  (6d) 


“  yil'^RaW 

where  m  is  an  integer  denoting  a  particular  dielectric  or  conducting  medium 
in  the  space  to  be  modeled.  We  also  define  the  proportiorial  electric-field 
vector 

E  =  R^E.  (7) 

Using  the  definitions  of  Equations  {6a)  -  (6e)  and  (7),  we  refO'*mulate  Yee's 

-V 

finlte-diffe'^ence  system  for  the  H  time-stepping  equations  as: 
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Assuming  m  “  MEDIA(i,J,k)  denotes  the  tvoe  of  medium  at  an  electric  field 


component  location,  _(1,j,k),  the  t  time-stepping  equations  are: 
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.  n+l  ^n  n+ij  n+^i 

Ej^d+ij.j.k)  »  Cg(m)E^d+ii.j.k)  +  Cjj(m)CH2d+>5.j-^,k)  -  H^d+^s. j-ij.k) 


n+Jj  n+^ 

+  Hyd+ij.j.k-Jj)  -  Hyd+!5,j,k+‘a)] 


(8d) 


n+l  n  n+‘s  n+'-j 

Eyd,j+!^,k)  =  C^(m)Eyd„j+‘j,k)  +  [H^d .  j+^,k+?i)  -  H^(i ,  j+!i,k-%) 


n+ij  n+H 

+  Hj,(i-i5,j+ia.k)  -  H^d+'i.j+'s.k)] 


(8e) 


,n+l 


n+*ii 


n+'n 


E^d.j.k+'i)  •  C^(m)E^d.j,k+hi)  +  C^(m)  [H^d+'i, j.k+^j)  -  H^{1-»a.j ,k+>i) 


n+^  n+^ 

+  H^(i  .j-^.k+^s)  -  H  (i[,j+*s,k+Hi)] 


(8f) 


ThiS  reformulation  eliminates  the  three  multiplications  needed  by  Yee 


to  compute  and  Further,  It  el Imlnates  the  need  for  computer  storage 


of  separate  e  and  a  arrays.  Now,  only  a  MEDIA  array  which  specifies  the 
type-integer  of  the  dielectric  or  conducting  medium  at  the  location  of  each 
electric  field  component  In  the  lattice  need  be  stored.  In  addition,  the 


and  a  of  each  medium  can  now  be  changed  without  having  to  re-punch  a 
large  data  card  deck,  if  the  basic  structure  geometry  is  unchanged.  Such 
a  change  involves  only  the  recalculation  of  the  few  values  of  Cg(m)  and 

Cjj(m). 

With  the  system  of  Equations  (8a)  -  (8f),  the  new  value  of  a  field 
vector  component  at  any  lattice  point  depends  only  on  its  previous  value 
and  on  the  previous  values  of  the  components  of  the  other  field  vector  at 
adjacent  points.  Therefore,  at  any  given  time  step,  the  computation  of  a 
field  vector  may  proceed  either  one  point  at  a  time,  or,  with  a  parallel 
processing  computer,  at  many  points  at  a  time. 

The  choice  of  iS  and  fit  is  motivated  by  the  reasons  of  accuracy  and 
algorithm  stability,  respectively.  To  Insure  the  accuracy  of  the  computed 
spatial  derivatives  of  the  electromagnetic  fields,  fi  must  be  Small  compared 
to  a  wavelength  (usually  <  A/10).  Further,  to  insure  that  the  cubic  lattice 
approximation  to  the  surfaces  of  the  structure  modeled  is  not  too  coarse, 

5  must  be  small  compared  to  the  overall  dimensions  of  the  structure. 

To  insure  the  stability  of  the  time-stepping  algorithm  of  Equations 
(8a)  -  (8f),  St  is  chosen  to  satisfy  the  Inequality 

■St  <  {-K  +  =  — - —  (for  a  cubic  lattice)  (9) 

"  6x^  5y^  6z^  c  /I 

uj,  ''max 

where  is  the  maximum  wave  phase  velocity  within  the  model  [9j.  The 
corresponding  stability  criterion  set  forth  by  Yee  In  Equations  (7)  and  (8) 
of  his  paper  [1]  is  incorrect. 

*‘^2,2,2  Lattice  Regions  and  Plane  Wave  Source  Condition  [14]-[16] 

As  shown  in  Figure  4a,  the  latest  formulation  of  the  FD-TD  lattice 
involves  the  division  of  the  computation  space  into  two  distinct  regions, 
separated  by  a  rectangular  surface  which  serves  to  connect  fields  in  each 
region.  In  two  dimensions,  the  surface  has  four  faces;  in  three  dimensions, 
the  surface  has  six  faces. 


16 


(a) 


Region  I  : 

Total  Fields  t 

io 

Region  2  : 
Scattered 
Fields 


•■’■■■P'"''-'® . 

i 


(b) 


Figure  4.  DIVISION  OF  FD-TD  LATTICE  INTO  TOTAL-FIELD 
AND  SCATTERED-FIELD  REGIONS,  (a)  Lattice 
division;  (b)  Field  component  geometry  at 
connecting  plane  y  “ 
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Region  I  of  the  FU-TD  lattice  is  denoted  as  the  total -f ield  region. 
Here,  it  is  assumed,  that  all  computed  field  quantities  are  comprised  of 
the  sum  of  the  incident  wave  and  the  scattered  field.  The  interacting 
structure  of  interest  is  embedded  within  this  region. 


Region  1  of  the  FD-TD  lattice  is  denoted  as  the  scattered-field  region. 
Here,  it  is  assumed  that  all  computed  field  quantities  are  comprised  only 
of  the  scattered  field.  The  outer  lattice  planes  bounding  Region  2,  called 
the  lattice  truncation  planes,  serve  to  implement  the  free-space  radiation 
condition,  discussed  in  the  next  subsection. 

The  rectangle  faces  comprising  the  boundary  between  Regions  1  and  2 
contain  E  and  H  field  components  which,  according  to  the  system  of  equations 
(8a)  -  (8f),  require  the  formulation  of  various  field-component  differences 
across  the  boundary  planes  for-  proceeding  one  time  step.  Typical  FD-TD 
computations  at  these  boundary  points  are  as  follows,  using  the  spatial 
coordinates  shown  In  Figure  4b. 


~n+1  ..n+1 

E  (i.j J  -  E  (i,J 


o^kqn.  (Hf)  ^ 


n+ij  n+V 


.•n 


(10a) 

(10b) 


-n+1 

Here,  usual  FD-TU  value  of  the  tp.tal_  component  evaluated 

at  point  (i,Jn)  and  time  step  n+1.  SimiVarlJf,  H  -ts)1s  the  usual 

FD-TD  value  of  the  scattered  H*  component  evaluated  at  point  (i.j^-^)  and 
time  step  ri+'j.  Cj^(tti)  denotes  a  proportionality  factor  defined  in  Equation 
(6e).  The  superscript  "i"  denotes  the  known  incident  field  component 
value.  These  computations  assure  consistency  of  the  subtraction  operations 
of  field  components  across  Llie  Region  1/Reglon  2  boundary.  In  effect, 
total-field  quantities  are  always  subtracted  from  similar  total -field 
quantltiesi  and  scattered-field  quantities  are  always  subtracted  from 
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similar  scattered-field  quantities.  This  enforcement  of  consistency  serves 
to  precisely  connect  the  total-field  and  scattered-field  regions.  Further, 
the  Inclusion  of  arbitrary  values  of  and  in  the  consistency  relations 
permits  the  precise  specification  of  the  desired  plane  wave  of  arbitrary 
angle  of  Incidence  and  arbitrary  polarization. 

There  are  a  number  of  key  advantages  to  this  methodology:  (a)  The 
high-dynamic-range,  total-field  formalism  Is  retained  for  the  entirety  of 
the  Interacting  structure,  permitting  accurate  computations  of  low-level 
fields  penetrating  Into  cavities  through  apertures,  and  In  the  shadow 
regions  of  scatterers.  (b)  The  scattered-field  formalism  Is  retained  for 
the  lattice  truncation  region,  permitting  a  very  accurate  simulation  of 
the  radiation  condition,  (c)  The  Incident  wave  contribution  need  be  com¬ 
puted  or  stored  only  for  the  field  components  at  the  rectangular  surface 
connecting  Regions  1  and  2.  This  results  In  much  less  computation  or 
storage  than  If  the  Incident  field  were  to  be  computed  at  all  points  within 
the  Interacting  structure  to  Implement  a  pure  scattered-field  formalism. 

(d)  The  scattered  near  field  in  Region  2  can  be  easily  integrated  to  derive 
the  far-fleld  scattering  and  radar  cross  section,  as  discussed  later. 

**2.2.3  Lattice  Truncation  Conditions 

A  basic  consideration  with  the  FD-TD  approach  to  solve  electromagnetic 
field  problems  Is  that  most  such  problems  are  usually  considered  to  be 
"open"  problems  where  the  domain  of  the  computed  field  Is  Ideally  unbounded. 
Clearly,  no  computer  can  store  an  unlimited  amount  of  data,  and  the  field 
computation  zone  must  be  limited.  The  computation  zone  must  be  large  enough 
to  enclose  the  structure  of  Interest,  and  a  suitable  boundary  condition  on 
the  outer  perimeter  of  the  computation  zone  must  be  used  to  simulate  the 
extension  of  the  computation  zone  to  Infinity.  Outer  boundary  conditions 
of  this  type  have  been  called  either  radiation  conditions,  absorbing 
boundary  conditions,  or  lattice  truncation  conditions  by  various  recent 
workers  in  this  area. 
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Let  US  now  consider  this  problem  In  more  detail.  Inspection  of 
Equations  (8a)  -  (8f)  Indicates  that  the  values  of  field  components  at 
lattice  truncation  planes  cannot  be  determined  from  the  system  of  finite- 
difference  equations  because  of  the  centered  nature  of  the  spatial  deri¬ 
vatives.  Therefoi'e,  these  values  must  be  computed  usinq  an  auxiliary  trun¬ 
cation  condition.  However,  care  must  be  taken  because  this  condition  must 
not  cause  the  spurious  back-reflection  of  outgoing  scattered  waves  (as 
observed  by  Yee),  and  must  not  cause  numerical  Instability.  The  goal  of 
formulating  the  truncation  condition  1s  to  make  the  lattice  truncation 
planes  Invisible  to  all  possible  waves  propagating  within  the  lattice,  as 
shown  In  Figure  1. 


A  desirable  truncation  condition  relates  In  a  simple  way  the  values  of 
the  field  components  at  the  truncation  planes  to  field  component  values  at 
points  one  oniiore  within  the  lattice.  We  now  review  several  possible 
approaches  which  have  been  recently  published. 


in  t^n^ee  dimensions,  an  outgoing 

scattered-wave  field  component,  F®  (either  an  electric  or  magnetic  field) 
has  a  (r,  0,i())  variation  of  the  type  [17]. 


(11) 


Here,  the  bracketed  Infinite  series  represents  In  effect  a  multi  pole  expansion 
of  the  scattered  field,  where  A,  B,  .  .  .  are  initially  unknown  functions  of 
angular  position. 

First-order  simulations  of  the  outer  lattice  boundary  condition  approxi¬ 
mate  the  A((),iii)/r  dependence  only.  These  approximations  include  those  of: 
a)  Taylor  et  al  [3],  using  a  radial  field  extrapolation  method;  b)  Taflove 
and  Brodwin  [9],  [12],  using  simulation  of  near-field  outgoing  waves  with 
propagation  delay,  a  field  averaging  process,  and  the  possible  use  of  elec¬ 
tric  or  magnetic  loss  in  the  scattered-field  zone;  and  c)  Merewether  [4]  and 
Kunz  and  Lee  [6]  simulating  the  radiation  condition  at  large  distances  from 
the  center  of  the  scatteror. 
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2-  Second  and  higher-order  approximations.  These  simulate  the 

O  ” 

B(0,(f)/r  and  higher-order  r-dependent  behavior  of  F®  In  addition  to  the 
A(o,i|>)/r  term.  These  approximations  Include  those  of:  a)  Engquist  and 
Majda  [18];  and  b)  Kriegsmann  and  Morawetz  [19].  In  general,  these 
approaches  involve  the  construction  of  mixed  radial  and  time-derivative 
operators,  on  F  ,  which  result  1n  D^(F  )  =  0  (r  )  so  that  arbitrary 
precision  at  any  r  can  be  obtained  by  Increasing  the  order,  1,  of  (here 
0  represents  the  order  of  a  function). 

A  second-order  approximation  of  this  type  at  truncation  plane  x  =  0 
Is  given  by  [18] 


L 

C  3x3t 


1 


.  „2pS  -2-5 

r  X  1  ^  X  9  r  \ 


X  »  0 


0. 


(12) 


for  the  radiation  condition  at  the  lattice  truncation  plane  x  =  0.  (Here 
c:  is  the  free-space  phase  velocity).  A  typical  FD-TD  computation  realizing 
this  radiation  condition  Is  given  by  [16]: 


-n+1 


(0,j,k+^J  =  -E5'’(l,j,k+%) 

t'2  (0,1+1  ,k+%)  -  2i5(0,j,k+Ji)  +  E;(0.j-l,k+^) 
+e"  (l.j+l,k+%)  -  2E"  (l,j,k+%)  +  E"(l.J-1,k+Js) 
+E;  (0,j,k+=^/t)  -2E;(0.j,k+%)  +  E"(0,J,k-'i) 


+e"  (l,j.k+Vz)  -2i;(1«J.k+%)  +  i;(l,j,k-Ji) 


(13) 


?>s'(wt+;sT 
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**n4*  1 

In  Equation  (13),  (0,j,k+%)  is  the  desired  truncation-plane  value  of  E 

at  point  (0,j,k+%);  iS  is  the  lattice  cell  size;  and  6t  is  the  time  step.  An 

analogous  condition  can  be  written  for  E^  at  x  =  0,  as  well  as  for  E-^  and 

E„  at  truncation  plane  ;  E  and  E  at  truncation  planes  y  =  0  and  y  ,  ; 
y  max  z  x  r  j  ^ 

and  E„  and  E„  at  truncation  planes  z  =  0  and  z  (in  three-dimensional 
X  y  ^  max 

problems). 

**2.2.4  Sinusoidal  Steady-State  Magnitude  and  Phase  Information 

Such  data  can  be  obtained  either  by  (a)  directly  programming  a  single¬ 
frequency  incident  plane  wave,  or  (b)  performing  a  separate  Fourier  trans¬ 
formation  step  on  the  pulse  waveform  response.  Both  methods  require  time¬ 
stepping  to  a  maximum  time  equal  to  several  wave  periods  at  the  desired 
frequency.  The  second  method  has  two  additional  requirements.  First,  a 
short-rise-time  pulse  suffers  from  accumulating  waveform  error 'due  to  over¬ 
shoot  and  ringing  as  it  propagates  through  the  space  lattice.  This  leads 
to  a  numerical  noise  component  which  must  be  filtered  before  Fourier  trans¬ 
formation.  Second,  Fourier  transformation  of  many  lattice-cell  field- 
versus-time  waveforms  would  significantly  add  to  the  total  requirements  for 
computer  storage  and  execution  time  [3,4]. 

Recent  IITRI  work  has  shown  that  very  accurate  magnitude  and  phase 
information  for  sinusoidal  steady-state  FD-TD  problems  can  be  obtained  by 
method  (a)  above  and  observing  the  peak  positive  and  negative-going  ex¬ 
cursions  of  the  fields  over  a  complete  cycle  of  the  incident  wave  (after 

having  time-stepped  through  2-5  cycles  of  the  transient  period  following 
the  beginning  of  time  stepping).  For  certain  two  and  three-dimensional 
scattering  problems,  a  dc  offset  of  particular  computed  field  components 
can  be  possible.  This  leads  to  the  following  requirements  to  obtain 
correct  magnitude  and  phase  data: 

1.  The  peak-to-peak  value  of  the  sinusoidal  response  at  any 
point  must  be  observed  to  eliminate  the  effects  of  any 
dc  offset  upon  the  computation  of  the  phasor  magnitude. 

2.  The  zero-crossing  of  the  field  waveform  may  not  be  useful 
in  determining  relative  phase.  Rather,  H  may  be  necessary 
to  locate  the  zero-derivative  points  of  the  waveform  for 
this  purpose. 
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The  following  is  an  algorithm  to  locate  the  zero-derivative  points  of 
a  field-vs-tirne  waveform  computed  using  the  FD-TD  method,  This  particular 
algorithm  has  an  uncertainty  of  one  time  step  or  less,  and  operates  by 
noting  where  the  sign  of  the  waveform  time  derivative  changes  from,  positive 
to  negative,  or  negative  to  positive. 

For  ease  of  understanding  the  algorithm  and  relating  it  to  the  Fortran 
program  of  Appendix  A,  the  algorithm  is  given  here  in  the  form  of  a  simpli¬ 
fied  Fortran  program.  The  symbols  used  here  are  defined  as  follows: 

NCYCS  =  the  number  of  complete  cycles  of  the  sinusoidal  Incident 
wave  that  the  FD-TD  program  is  time-stepped; 

FREQ  =  the  frequency  of  the  sinusoidal  incident  wave  (in  Hz); 

DT  =  the  time  step  of  the  FD-TD  program  (in  seconds); 

NMIN  "  the  number  of  time  steps  spanning  one  complete  cycle  of 

the  sinusoidal  incident  wave; 

NMAX  “  the  total  number  of  time  steps  that  the  FD-TD  program  is 
time-stepped; 

NMINA  »  the  number  of  time  steps  spanning  one-half  cycle  of  the 
sinusoidal  Incident  wave; 

N  »  absolute  time-step  number  (from  beginning  of  time-stepping), 
ranging  from  1  to  NMAX; 


NA  »  number  of  time  steps  since  the  start  of  the  present  full- 
wave-cycle  observation  period  for  zero-derivative  points, 
ranging  from  0  to  NMIN-1 ; 

F  =  field  component  (either  E^,  E^,  E^,  or  H^); 

DK,  DFNU  »  time  derivatives  of  the  F  waveform; 

ENF  =  stored  value  of  F,  taken  at  the  time  of  each  zero- 
derivative  event  of  F; 

IPHF  =  stored  value  of  NA,  taken  at  the  time  of  a  zero-derivative 
event  of  F  coupled  with  a  change  of  sign  of  the  F  derivative 
from  positive  to  negative. 

The  algorithm  is  as  follows: 
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NMIN  =  I .U/hRbg/UI 

NMID  «  NMIN . 

NMAX  =  NCYCS  *  NMIN 

NMINA  “  FL0AT(NMIN)/2.0 

NMIDA  «  NMINA 

NMAXA  =■  NCYCS  *  2  *  NMINA 

DO,  21  ,  NMAX 

FNMOD  =  FLOAT (N)/FLOAT(NMID) 

NA  *  N  -  IFIX(FNMOD)*NMID 
DO  Z22  J  =  .lUmlts 
DO  zzz  I  -  iHiiiits 
STORE  F(I,J) 

F(I,J)  '  new  F  computed  using  Eqns.(8a-8f) 
DFNU  =  F(I,J)  --  STORE 
IF(OFNU-*DrU,J).6T.O.O)  GO  TO  xxx 
ENF(I,J)  «  Fd.J) 

IF  (DFNU. LT. 0.0)  IPHF{I,0)«  NA 
xxx  DF(I,J)  =  UFNU 
ZZZ  CONTINUE 


Define  constafits 


Test  field  component 
F(1  J)  for  change  of 
sign  of  time  derivative 


no  12  K^^NMINA, NMAXA, NMIDA 
IF(N.EQ.K)  GO  TO  13 

12  CONTINUE 
GO  TO  20 

13  Fortran  statements  which  first  print  out 
the  ENF  array  and  then  reset  the  ENF  array 
to  zero. . 


Test  for  half -cycle 
Intervals  of  time- 
stepping 


20  DO  512  K*NMIN.NMAX,NM1D 
IF  (N.EQ.K)  GO  TO  513 
512  CONTINUE 
GO  TO  520 

51 J  Fortran  statements  which  first  print  out 

the  IPHF  array,  and  then  reset  the  IPHP  array 
to  zero. 


Test  for  full -cycle 
Intervals  of  time- 
stepping 


b«iU  CUNllNUL 
21  CONTINUE 
STOP 
END 

Every  tOne  step,  the  above  algorithm  tests  the  desired  field  waveform, 
P,  for  a  change  of  sign  of  the  time  derivative.  If  a  sign  change  is  de¬ 
tected,  the  latest  computed  value  of  F  is  stored  in  the  array  ENF.  Note 
that  this  value  can  be  either  positive  or  negative.  If,  further,  the  new 
sign  of  the  time  derivative  Is  negative,  indicating  that  a  "peak"  of  the 
F  waveform  has  just  occurred,  the  tinie‘-step  number,  NA,  Is  stored  1n  the 
array  IPHF. 

The  above  algorithm  also  tests  to  determine  when  half-cycle  and  full- 
cycle  Intervals  of  time  stepping  have  passed.  Each  half-cycle,  the  stored 
values  of  ENF  are  first  printed  out  and  then  reset  to  zero.  Each  full- 
cycle,  the  stored  values  of  IPHF  are  first  printed  out  and  then  reset  to 
zero.  Therefore,  the  desired  zero-to-peak  value  of  F  Is  determined  by 


End  of  main  time- 
stepping  loop 


subtracting  the  last  printed  value  of  ENF  from  the  immediately-preceding 
printed  value  of  ENF,  taking  the  absolute  value,  and  then  dividing  by  two. 
Also,  the  desired  relative  phase  of  F  can  be  obtained  by  subtracting 
values  of  IPHF  (as  given  in  the  last  printout  of  IPHF)  from  some  fixed 
‘reference  value  at  a  specified  point.  This  methodology  avoids  any  ambiguity 
due  to  either  a  possible  DC  offset  of  the  F  waveform  or  the  repetitive 
nature  of  the  F  waveform. 

**2.2.5  Far-Fleld  Scattering  Information  via  the  Near-to-Far 
Field  Transformation 

In  principle,  the  electromagnetic  scattering  by  an  arbitrary  conductiiu'i 
body  can  be  determined  by  solving  an  integral  equation  for  the  induced 
electric  currents  on  the  surface  of  the  body.  Then,  the  induced  current^ 
can  be  integrated  to  calculate  the  near  or  the  far  fields.  However,  che  ' 
body  surface  may  have  a  complex  shape  or  may  be  loaded  with  dielectrics  in 
some  way  as  to  make  implementation  of  the  needed  integral  equation  and 
surface  integrals  very  diffjcuj  t,  and  in  fact,  a  unique  problem  for  each 
scafterer.  A  useful  alternative  would  be  to  obtain  the  scattered-field 


_ _  * 
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information  from  off-surface,  near  field  data,  rather  than  surface-current 
data.  Here,  the  near  field  information  could  be  computed  using  the  FD-TD 
method,  which  can  easily  account  for  the  complexities  introduced  by  the 
object's  shape  and  composition.  Further,  the  near-field  data  could  be  inte¬ 
grated  along  arbitrary,  planar,  virtual  surfaces  which  completely  surround 
the  object  of  interest.  In  this  manner,  accounting  of  the  data  would  be  simple 
and  would  be  independent  of  the  precise  shape  of  the  object  which  resides 
within  the  integration  surfaces.  Fortunately,  this  near-field  to  far-field 
transformation  can  be  shown  to  be  feasible  by  using  powerful  electromagnetic 
equivalence  relations. 

f 

Figure  5  demonstrates  a  procedure  to  implement  the  electromagnetic 
equivalence  needed  to  obtain  far-field  data  from  FD-TD-computed  near-field 
information  at  an  arbitrary,  closed,  virtual  boundary,  S  ,  fully  enclosing  a 

d 

scatterer.  S,  separates  Region  A  (FD-TD  computation  zone)  from  Region  B  (ex- 

a 

terior  scattered-field  zone).  is  optimally  a  rectangular  surface  which 
conforms  with  distinct  planes  of  the  FD-TD  lattice  located  in  the  scattered- 
field  region  (Region  2)  of  Figure  4. 

To  apply  the  hybrid  FD-TD  method,  the  tangential  components  of  the  scat¬ 
tered  fields,  and  are  first  obtained  at  S  using  FD-TD.  Then,  as  indi- 

d 

cated  in  Figure  5b,  an  equivalent  problem  is  set  up  external  to  S  which  is 

ct 

completely  valid  for  Region  B.  The  new  excitation  data  are  J  and  ,  the 

^eq  ""eq 

equivalent  surface  electric  and  magnetic  currents,  respectively,  on  S  which 

d 

are  obtained  according  to  [20]; 

^  (f)  =  n  X  ff^(r)  (14a) 

eq 

(r)  =  -n  x  P{r)  (14b) 

eq 

A. 

where  n  is  the  outward  unit  normal  vector  at  the  surface  S  . 

d 

The  equivalent  surface  currents  on  produce  the  same  scattered  field 
(f^,  it^)  external  to  as  in  the  original  problem.  Region  A  is  now  made 
empty  with  zero  fields  and  no  sources,  and  also  forced  to  have  the  same 
medium  characteristics  as  Region  B. 
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Figure  5.  NEAR-TO-FAR  FIELD  TRANSFORMATION  GEOMETRY 
(d)  Original  problem;  (b)  Equivalent 
problem  external  to  the  virtual  surface,  S. 

Ql 
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The  scattered  far  fields  are  thus  given  by  the  transform  of  the  equiva- 
lent  currents  of  Eqns.  (14)  over  the  free-space  Green's  function  [17],  If 
IPq.  Eqi  Oq)  are  the  Region-B  medium  characteristics,  and  if  we  define  the 
relation  s  *  jui,  we  have  the  following  general  scattered-field  expressions 
[17]  in  three  dimensions  for  r  >  r  : 

d 


^  t  Viv/(^(r\s)l  -  •V^„(r.s}  }  + 

Vx^(r.s) 


(15a) 


"(^(r.s)  -  _s_  {  V[V-/((r,s)]  "Y„^A(r,s)}  + 

2  ° 

^0 

(15b) 

T5'~+"si;i 

In  Eqns.  (15a)  and  (15b),  the  electric  vector  potential  is  given  by 


'“o  *  “"o’ 

4  F  s 


K  G{r,r/,s)  dS/ 

*eq  *  ®  * 


(16a) 


G(r^r^s) 


-Yjr-r;' 


I  “>• 

r  -  r* 


(16b) 


I 


[  (x-x/)"  +  (y-y/)^  +  (t-z /)^ 


..»2  ,>ii 


(16c) 


and  the  magnetic  vector  potential  is  given  by 
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(16d) 


(16e) 


Taking  the  limit  r  ■>  «  gives  the  required  far-fleld  scattering  distribution 
1171. 

in  order  to  validate  the  feasibility  of  the  proposed  hybrid  method,  the 
scattering  of  a  transverse-magnetic  (TM)  polarized  plane  wave  by  a  two-dImen 
sionai  conducting  cylinder  of  arbitrary  cross  section  was  considered  during 
the  present  research  program.  This  canonical  problem,  shown  in  Figure  6,  Is 
well  documented  [17,  21,  22]  and  therefore  serves  as  a  good  test  of  feasibll 
ity  and  accuracy  of  the  hybrid  method. 

The  near-field  to  far-fleld  transformation  discussed  above  Is  now  spec¬ 
ialized  for  this  two-dimensional  canonical  problem.  Assuming  that  Is  the 
free-space  propagation  constant  and  Is  the  Intrinsic  Impedance  equal  to 
we  have  In  the  far-fleld  region: 

‘eq  ^eq  '^eq 


where 


■j  ^  ,  k  K  f  I-'-  ill-  (17b) 

^eq  0  0  J  «q. 


r  ^  k  ,  f  («'.  y')  *  l-sln*) 

”  J  '"Ik.y 


K  -  § _  e"J  T 

/Bmkp 


29 


CROSS  SECTION 


Figure  6.  GtOMtTRY  OF  TWa-DlMFi'iSIONAL  FCATFERING,  TM  POLARIZATION  CASE 
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and  the  scattering  cross  section  based  on  the  equivalent  currents  is  given 
by: 

RCS„„  «  11m  {  Zvp  |E®(<t))/EV  }  •  (17e) 

P  oo 

For  comparison,  the  far  fields  and  scattering  cross  section  obtained  by 
the  usual  surface  electric  currents  (J^,)  approach  Is  based  upon  solving  the 
following  Integral  equation  using  the  method  of  moments  (MOM): 


J2(P^)  dr 


for  p  E  C 


(18a) 


eJ(p) 


I18b) 


In  Eqns.  (Ifla)  and  (10b),  Is  the  angle  of  Incidence  (angle  between  the  x 
axis  and  the  dircctlort  of  propagation),  and  Is  the  Hankel  function  of 
the  second  kind  and  zero  order.  Now,  the  far-fleld  distribution  Is  obtained 
directly  from  the  Induced  cylinder  currents,  0^,  by: 


^  l:2\.p-J(kP  +  ^^)  f  J,(X',  y^) 

C 


(19) 


The  scattering  cross  section  Is  again  given  by  Equation  (17e). 

2.2.6  Penetrating  Near-F1eld  Information  via  the  Schelkunoff 
Aperture  Electric  Current  Equivalence  Principle  [13] 

The  analysis  of  the  electromagnetic  excitation  of  an  aperture  on  an  arb¬ 
itrarily  shaped  object  Is  a  complex  problem  [231.  This  problem  has  been  given 
special  attention  In  the  broad  area  of  EMP  penetration  and  simulation  studies 
[241  with  great  efforts  being  applied  to  assess  coupling  to  objects  present 
behind  apertures  on  finite,  metallic,  hollow  scattering  bodies. 


The  problem  of  the  penetration  and  coupling  of  electromagnetic  energy 
through  an  aperture  has  been  studied  extensively  by  many  investigators  using 
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analytical  and  iiietiio't  I'l  •■muiiionts  (MOM)  numet'icdl  approachos  l;M|  to  charac¬ 
terize  the  behavior  of  simple  aperture  shapes  In  a  conducting  screen  [25]  - 
[29]  or  in  a  finite,  scattering  body  [30,  311.  However,  the  analysis  becomes 
very  complex  1f  there  are  other  scattering  objects  In  the  vicinity  which  are 
coupled  to  the  aperture  [23],  [32,  33]  wherein  the  mutual  Interaction  must  be 
fully  taken  into  account. 

One  powerful  alternate  approach  Is  the  pure  FD-TD  method  previously  dis¬ 
cussed  In  Re  '  . ,  trices  [12]  and  [13]  which  allows  computation  of  the  penetrated 
Internal  electromagnetic  fields  by  direct  modeling  of  the  complete  object  of 
Interest,  Including  apertures  and  Internal  contents.  However,  the  pure  PD-TD 
method  Is  best  suited  for  modeling  localized  regions.  When  such  a  region  Is 
part  of  a  larger  structure  that  significantly  affects  the  electromagnetic  ex¬ 
citation,  It  becomes  difficult  to  account  for  the  physics  of  the  entire  coup¬ 
ling  problem  using  a  single  FD-TD  model  having  a  constant  space  lattice  cell 
size.  In  fact,  electromagnetic  coupling  problems  Involving  the  need  to  account 
simultaneously  for  the  effects  of  both  large  and  small  spatial  details  (the 
"global -local  problem")  have  been  difficult  to  structure  with  any  single  ana¬ 
lytical  or  nunierical  technique. 

In  order  to  treat  such  coupling  problems  more  effectively,  a  new  hybrid 
MOM/FD-TD  technique  was  developed  based  on  a  novel  use  of  a  field  equivalence 
theorem  due  to  Schelkunoff  [13],  [20],  [34  ,  35],  This  hybrid  technique,  des¬ 
cribed  in  this  section,  is  basically  an  equivalent  aperture  excitation  method. 
This  allows  analysis  of  the  coupling  problem  In  two  distinct  steps: 

Step  1-  Analysis  of  the  relatively  simple  exterior  problem  using 
MOM  to  compute  the  equivalent  excitation  currents  in  the 
apertures  leading  to  the  Interior  region.  This  can  be 
done  independent  of  any  knowledge  of  the  contents  of  the 
interior  region . 

Step  2"  Use  of  the  FD-TD  method  to  analyze  the  relatively  complex 
interior  region,  assuming  as  an  excitation  the  equivalent 
currents  found  In  Step  1. 

In  this  way,  each  analysis  method  is  applied  1n  the  range  of  structure  size 
and  complexity  that  it  Is  best  suited  for,  allowing  an  overall  solution  that 
Is  accurate  for  large,  simple  structures  that  are  penetrated  by  apertures 
leading  to  complex,  interior  cavities. 
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Figure  7a  illustrates  the  classical  problem  of  a  perfectly  conducting 
scatterer,  ,  with  an  aperture,  S^.  External  to  the  scatterer  in  Region  1, 
an  incident  field  illuminates  the  aperture,  and  part  of  the  energy 

penetrates  into  the  cavity  (Region  2)  through  the  aperture.  To  compute  the 
total  field  Region  2,  many  analytical  and  MOM  numerical  approach¬ 

es  studied  in  the  literature  have  employed  integral-equation  formulations 
treating  Regions  1  and  2  simultaneously.  These  approaches  characterize  either 
or  in  terms  of  tangential  fields  or  equivalent  surface  currents  [20], 
[211 . 

An  alternative  is  to  employ  a  field  equivalence  theorem  due  to  Schelkun- 
off  (Theorem  no.  3  in  [20])  to  permit  sequential  treatment  of  Region  1  and 
Region  2.  This  theorem  is  based  upon  an  equivalent  aperture  electric  current 
formulation  which  connects  the  exterior  and  interior  regions,  This  formula¬ 
tion  expresses  fields  in  a  region  as  the  superposition  of  the  so-called 
short-circuit  fields  (^gg*  with  the  aperture  not  present  (shorted)  plus 
the  aperture  field  contribution,  (f^,  rf^),  maintaining  the  required  continui¬ 
ty  of  the  fields  across  the  aperture.  The  first  partial  field,  i^,.),  is 

SC  S  w 

simply  equal  to  zero  in  the  Region  2.  The  second  partial  field,  designated 
(^^,  ll^),  is  generated  by  the  non-physical  current  distribution,  flowing 
through  the  empty  space  of  the  open-circuit  aperture  locus,  where 
the  short-circuit  aperture  current  distribution. 

In  Figures  7b  and  7c,  the  MOM/FD-TD  hybrid  method  is  illustrated  sche¬ 
matically,  and  is  summarized  below  as  a  four-step  computation. 

Region  1:  MOM  Technique 

(a)  The  aperture  region,  S^,  is  short-circuited,  and  the  straight¬ 
forward  exterior  problem  is  solved  via  MOM  to  obtain  the  in¬ 
duced  electric  current  distribution,  in  the  short-circuited 
aperture  region  (Figure  7b). 

(b)  is  now  placed  in  the  open  aperture  region  with  a  sign  change 
(Figure  7c)  to  account  for  the  continuity  of  the  fields  in  the 
aperture  region.  The  electric  current,  non-physically 
placed  in  the  open  aperture  region  acts  as  an  equivalent  source 
for  the  fields  in  Regions  1  and  2. 
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REGION  ©(f|,fii) 

A 


MOM 


Figure  7.  HYBRID  MOM/m-Tl)  TECHNIQUE  FOR  COMPUTING  PENETRATION 
INTO  CAVITIES,  (u)  Original  problenii  (b)  MOM  solution 
of  exterior  problem  for  shorted-aperture  case*,  (c) 
FD-'TD  solution  of  Interior  problem  using  equivalent 
aperture-cu Trent  excitation 
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Region  2:  FD-TD  Technique 

(c)  We  note  that  the  MOM  technique  gives  the  equivalent  current 
source  distribution,  In  the  frequency  domain  as  a  spatial 
distribution  of  phasor  quantities  having  both  magnitude 
(relative  to  the  incident  fields)  and  phase  (relative  to  some 
phase  reference,  normally  at  the  origin  of  the  coordinate 
i.ystem).  With  the  FD-TD  method  being  a  time-domain  technique, 
the  phase  distribution  of  -5^,  Is  Interpreted  as  a  time-delay 
distribution  with  respect  to  the  original  phase  reference  lo¬ 
cation.  The  magnitude  distribution  Is  taken  Intact  without 
interpretation  or  modification.  In  this  manner,  the  FD-TD 
aperture  equivalent  current  source  distribution  assumes  sinu¬ 
soidal  steady-state  quantities  starting  at  the  very  beginning 
cf  time-stepping,  with  the  proper  time  delay  to  account  for 
phase  shift. 

(d)  Using  the  FD-TD  approach,  the  Interior  fields  (?2,  1^2^ 
puted  directly  by  using  as  a  source  term  distribution  In 
the  V  X  H  difference  equations.  For  example,  If  -0^  Is  a 
source  term,  Eqn.  (8d)  Is  re-written  for  the  apertfire  points 
as; 

m  =  MEDIA  (1+Js,  J,k) 
i"^^(1+^.j,k)  »  Cg(m)ij(l+is,j,k)  + 

-H;^‘*(1H,,:P),,k)  +  (1+h;.j,k-»5)  -  Hy  {i+h,iMh) 

n+hs 

-Jq  (1+»a.j.k)]  (20) 

X 

Here,  specifying  in  the  aperture  Is  equivalent  to  sped-  • 
fylng  the  addition  ’^of  a  d1scont;1nu1ty  in  the  z-directed 
tangential  magnetic  field,  H^,  across  the  aperture  source  plane. 
Namely,  we  have  added 
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where  and  are  tangential,  magnet iu  (ieldi.  located  an 
itifinitesimaV  distance ‘to  either  side  of  the  equivalent 
aperture  source  plane.  This  procedure  is  consistent  with 
the  partial -field  approach  discussed  by  Schelkunoff  [20]. 

Note  that,  unlike  possible  hybrid  MOM/MOM  methods,,  the  hybrid  MOM/ 

FD-TD  method  does'not  require  computation  of  an  equivalent  aperture  electric 
field  excitation,  1-^.  This  is  because  the  interior-region  FD-TD  solution 
easily  accepts  the  non-physical  aperture  electric  current  distribution, 

-Jp,  as  the  excitation.  Thus,  there  is  no  need  to  set  up  and  solve  for 
the  mutual  interactions  of  the  cavity  contents  and  the  apertures,  and  no 
need  to  compute  the  cavity  Green's  function.  The  hybrid  MOM/FD-TD  method 
easily  permits  the  modeling  of  the  cavity  interior  to  as  fine  an  extent  as' 
desired  without  any  additional  numerical  complication.  For  realistic, 
general  cavities  having  Green's  functions  that  are  difficult  or  Impossible 
to  compute,  the  hybrid  MOM/FD-TD  approach  may  be  the  only  way  to  calculate 
the  penetrated  interior  fields. 

2.2.7  ThJ_n -Wi  re  Co_upJJ 

2 . 2 . 7 . 1  Simplified  Model  for  Total -Field  FD-TD  Formulation 

A  very  simple  way  of  approximating  electromagnetic  coupling  tn  a 
round,  thin  wire  having  a  diameter,  2r  ,  less  than  the  FD-TD  lattice-cell 
size,  iS,  Is  Indicated  in  Reference  [13].  This  method  can  be  easily  used 
when  the  wire  of  interest  is  embedded  in  a  total -field  modeling  region  of 
the  FD-TD  lattice.  It  1s  based  upon  modeling  the  thin  wire  as  a  thicker, 
square-cross  section  wire  occupying  exactly  one  cell  of  the  lattice.  Simply, 
the  electrical  conductivity,  o^,  assigned  to  the  E  components  tangential  to 
the  model -wire  space  cells  is  adjusted  so  that  the  high-frequency  internal 
impedance  of  the  model  wire  equals  that  of  the  actual  thin  wire  of  interest. 
For  example,  the  high-frequency  Internal  Impedance  of  a  thin  wire  is  given 
by  [36] 
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‘  1 


Wire 


2'nr , 


ohms /meter 


(22a) 


where  f  Is  the  frequency  of  the  wave  excitation,  Is  the  permeabnity 
of  the  wire  material,  1s  the  conductivity  of  the  wire  material,  and  r 
is  the  radius  of  the  wire.  The  high-frequency  Internal  impedance  of  a 
model  wire  having  a  square  cross  section  occupying  one  lattice  cell  Is 
well  approximated  by 


w 


flTfU, 


model  ^  0,,, 
wire  ' 


(1+i)  ohms/ meter 
''46 


(22b) 


where  Is  the  vacuum  permeability,  Is  the  conductivity  of  the  model 
wire  material,  and  iS  Is  the  width  of  the  model  wire  (one  lattice  cell  size), 
Equating  {22a)  and  (22b),  the  desired  conductivity  of  the  model  wire  Is 
given  by 


(irr^/26) 


nihos/m 


(22c) 


It  should  be  pointed  out  that  this  procedure  achieves  extreme  sim¬ 
plicity  at  the  expense  of  neglecting  the  exact  geometrical  relation  of  the 
thin  wire  and  Its  surroundings.  The  Internal  impedance  of  the  thin  wire  Is 
adequately  modeled  here,  but  the  precise  capacitance  and  inductance  linking 
the  thin  wire  to  adjacent  structures  is  somewhat  In  error  due  the  larger 
effective  radius  used  for  the  model  wire.  The  following  discussion  sum¬ 
marizes  a  procedure  for  taking  Into  account  the  radius  of  the  thin  wire 
with  regard  to  these  external  Impedances.  This  procedure  Is  considerably 
more  complicated  than  the  above,  and  is  usable  mainly  for  a  scattered- 
fleld  version  of  the  FD-TD  method.  The  relative  accuracy  of  the  two 
approaches  remains  to  be  determined, 
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**;^.2.7.2  Holland-Simpsoh  Model  [8]  for*  Scattered-Field 
FD-TD  Forynglatlon  _  •  ; '  ^ 

Holland  and  Simpson  [8]  have  recently  described  a  scattered-field- 
only  FD-TD  code  (which  they  name  "THREDE")  that  permits,  the  modeling  of 
electromagnetic  coupling  to  thin  wires  having  a  diameter  of  less  than  one 
space  cell  of  the  FD-TD  lattice.  This  sub-section  summarizes  the  theoreti¬ 
cal  background  and  the  tlrln-wire  couoling  algorithm  published  by  these 
authors. 


Figure  8  shows  a  cross  section  of  a  FD-TD  cell  with  a  wire  running 
perpendicularly  to  the  plane  of  the  figure.  Holland  and  Simpson  separate 
the  electroinagnetic  field  into  an  incident  and  a  scattered  part 

t  -  i  (23a) 


•  >c;  •  ^  ^ 

H  "  H^  +  H’ 


(23b) 


They  then  write  Maxwell's  equations  in  polar  coordinates  about  the  wire 
center,  ’’'he  0  component  of  the  curl-E  equation  is  given  by 


PE  PE  Of  'at  ■' 
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hz  uz  ;ir  'dr 
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At  r  a,  there  is  the  boundary  condition 


s  i 


They  integrate  Equation  (24)  out  to  some  radius  r  to  yield 


E,  (r)  +  E,'(r)  . 


t  I  (t,^  +  E,')dn 

a 


Figure  8.  GEOMETRY  FOR  COMPUTING  THE  IN-CELL  INDUCTANCE 

OF  A  WIRE  IN  A  CELL  OF  RECTANGULAR  CROSS  SECTION  [3] 
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Figure  9.  PARTITIONING  I'’  INTO  CURRENT  DENSITIES  AT  THE 
FOUR  .CLOSEST  E®  MESH  POINTS  [8] 
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By  hypothesis,  the  cell  sizes,  and  thus  all  r's  of  Interest,  are 
small  compared  to  the  shortest  wavelength  present.  Then,  in  an  average 
sense  around 0,  Holland  and  Simpson  write 


(27) 


r 


.  . 

2nrt„ 


(28) 


where  I  is  the  current  on  the  wire  and  Q  is  the  charge  per  unit  length  on 
the  wire.  If  the  wire  is  driven  by  current  injection,  locally  represents 
the  relevant  injected  current  (Normally,  a  scatterer  would  not  be  excited 
by  an  incident  field  and  current  injection  inthe  same  event.  Thus,  Holland 
and  Simpson  usually  have  either  E"*  0  or  '0).  The  and  Q®  represent 

the  scattered  current  and  charge  per  unit  length  on  the  wire.  Holland  and 
Simpson  use  Equations  (27)  and  (28)  to  re-express  Equation  (26)  as 


E/‘(r)+E^^r) 


0(I®+I^] 

■  li't . 


+  31.1  •  JnjrZal 
az  SttEq 


(29) 


They  provide  similar  treatment  of  the  r  component  of  the  curl  H  equation, 
resultint)  in 


<:  c  H 

ag-’  .  ;.)i  ..  m 
at  nz  "  nz 


(30) 


(Integrating  over  u  fi’oiii  0  to  2ir  to  obtain  the  average  around  0  causes  the 
(l/r)(i)H^/iiii)  term  in  the  cur1-H  equation  to  go  to  zero.  Equation  (30) 
alternatively  can  he  interpreted  as  merely  a  statement  of  conservation  of 

charge  on  the  wi re. ) 
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Holland  and  Simpson  point  out  that  the  Incident  quantities 
h\  and  are  known,  externally-defined  variables  and  are  input  quantities 
to  their  THREDE  FD-TD  system.  However,  they  note  that  the  scattered 
quantities  E^,  Q®,  and  I®  are  unknowns  which  are  solved  for  by  means 

of  their  FD-TD  program. 

They  next  re-write  Equation  (29)  as 
E/(,r)+E,^(r) 

L  1  Z 


0,  r  <  a 


They  state  that  most  of  the  uncertainty  in  their  modeling  of  the  coupling  to 
a  thin  wire  depends  on  the  next  operation  performed  on  Equation  (31).  Their 
common  manipulation  Is  to  average  Equation  (31)  over  the  area  bounded  by 

p 

r  «  R  where  itR  •  the  cell  area,  AX^'AY^.  Averaging  removes  the  r  dependence 
of  the  equation  and  yields 


L 


nt 


S'^o 


*  <E 


(32) 


where  L  has  the  dimensions  of  Inductance  per  unit  length  and  is  given  by 


i 


S 

i 

I 


j 

I 
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s. 

g. 

T 
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4 


Hollantl  and  binipsoti  found  eiiipiriLa'Ily  [8]  Lhat  the  lower  limit  of 

r=0  on  the  denominator  Integral  of  Equation  (33)  results  in  a  better  fit 
of  their  computed  results  to  known  analytic  solutions  of  simple  problems 
than  does  a  choice  of  r=a.  Thus,  the  limit  choice  used  1n  Equation  (33) 
is  more  pragmatic  than  rigorous.  The  numerator  integrisl  of  this  equation 
must  begin  at  r=a  because  vanishes  for  r  <.  a  as  Equation  (31) 

indicates. 

Holland  and  Simpson  implement  their  thin-wire  coupling  model  by 
applying  a  finite-difference  analog  of  Equations  (30)  and  (3i?)  as  well  as 
Maxwell's  equations.  Thus,  a  cell  with  a  wire  running  through  it  requires 
eight  quantities  (six  field  comtJOi.u.,ts  and  and  I®)  to  be  advanced  each 
program  time  step.  For  a  wire  segment  parallel  to  the  7  axis,  they  found 
that  it  is  easiest  computationally  to  place  mesh  points  in  the  planes 
of  mesh  points,  aii-i  mesh  points  in  the  planes  of  mesh  points. 

The  Hoi  1  and-Simpson  th'"  wire  coupling  algorithm  thern'fore  proceeds 
as  follows  L8]: 

1)  Evaluate'.!^  ><md  “•  l'  at  the  I"*  mesh ‘points. 

2)  Evaluate'  ■•at,tiie  mesh  points.  For  a  wire  along  z,  It 

is  desirable  to  run  the  wire  through  the  mesu  points.  Then 
<  E^^ >  is  just  E^'  at  that  point.  If  the  wire  cannot  conveniently 

be  run  through  the  E^^  mesh  point,  <E  should  be  Interpolated 

^  * 

from  the  four  closest  mesh  points,  as  shown  1n  Figure  t;  or  9. 

b)  Advance  according  to  the  following  finite-difference  analog  of 

Equation  (32); 


(34) 
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4)  Advance  in  the  usual  way  from  the  curl-T  equation. 

5)  Advance  time  by  At/?. 

6)  Partition  I®  +  among  the  four  closest  mesh  points  according 

to  Figure  9.  If  the  wire  runs  throug'ti  the  E_®  mesh  point,  as  pre- 
si  *  s 

ferred,  all  I  +  r  are  associated  with  that  mesh  point,  {I,J,K)» 

(I®(K)  +  1‘'(K))/(aX^.AY^). 

7)  Advance  E^  using  the  curl-H  equation  with  0®  as  determined  from 
partloiiing  the  I®  +  according  to  step  6. 

8)  Advance  0®  according  to  the  following  finite-difference  analog  of 
Equation  (30): 


q®(K)'^^^  *  Q®(K)^"^ 


[®(Kl"+l^(Kl"-I®(K-ll'’-I^(K-ll" 


9)  Advance  time  by  At/2  and  go  back  to  step  1. 

Holland  and  Simpson  [8]  also  describe  their  treatment  of  wire  Junctions; 
evaluation  of  the  In-cell  Inductance  for  non-cubic  cells;  boundary  conditions 
for  wires  ending  either  In  free  space  or  at  a  conducting  surface;  resistive 
wires;  and  applications  to  a  linear  antenna  and  a  circular  loop  antenna.  The 
reader  1s  referred  to  their  paper  for  the  details. 


3.0  EXAMPLES  OF  COMPUTED  RESULTS  OF  THE  FD-TD  METHOD 

3. 1  Pure  FD-TD  Method,  Three-Dimensional  Penetration  Problems 

3.1.1  Empty  Cylindrical  Metal  Cavity,  Broadside  Incidence, 

Transverse  Electric  (TE)  Polarization  Case  [13].  _ 

This  problem  Involved  the  computation  of  the  fields  within  a  19.0-cm 
diameter,  68.5-cm  long,  circular  aluminum  cylinder  with  one  open  end.  The 
incident  plane  wave  was  assumed  to  have  a  frequency  of  300  MHr.,  have  the 
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field  compoiitiiiLi  L  and  li  ,  and  propayat.”  perpc'ndicul ar  tu  tho 
■  inc  ^Inc 

cylinder  axis  with  a  transverse  electric  (TE)  polarisation.  Since  the 
cutoff  frequoncy  of  the  cylinder  (as  a  waveguide)  Is  above  900  MHz,  the 
Interior  fields  are  expected  to  decay  with  distance  from  the  aperture. 

The  geometry  of  the  cylinder  inode!  relative  to  the  160  x  63  x  24-0611 

(1,450,000  unknown  field  components)  FD-TD  problem  lattice  is  depicted  in 

Figure  10  at  lattice  pynmetry  plane  k  -  24  and  in  Figure  11  at  a  typical 

1  *  constant  plane  cutting  through  the  cylinder  transversely.  Like  the 

ax1dl“1nc1denc6  case  discussed  In  [11],  4  *  0.5  cm  «  Xjj/200  was  used; 

6t  “  8.33  psec;  and  aluminum  a  *  3.7  '  1o’ mho/m.  First-order  correct 

radiation  conditions  were  available  at  the  time;  thus,  an  air  conductivity, 

a  --  0.01  mho/m,  was  used  to  attenuate  scattered  waves  and  accelerate  con- 

£ 

vergence  of  the  Interior  fields.  2.5  ‘  10  words  of  memoy'y  and  9.4  minutes 
of  central  processor  time  were  required  on  the  Control  Data  Star  100  for 
a  complete  run  of  800  time  steps  (2.0  periods  of  the  Incident  wave  at 
300  MHz).^ 

Figure  12  compares  the  FD-TD  computed  ]  along  the  cylinder 

^1nc 

axis  with  computed  results  by  D.  Wilton  and  A.  Gllsson  using  a  frequency- 
domain,  MOM,  body-of“revolut1on  computer  program  [37],  The  FD-TD  results, 
shown  as  a  solid  curve,  are  after  800  time  steps;  the  MOM  results  are 
graphed  as  a  dashed  curve.  Excellent  agreement  (within  0.5  dB)  of  the 
results  of  the  two  approaches  1s  seen  over  the  first  15  cm  of  field  pene¬ 
tration  Into  the  cylinder.  Over  this  span,  the  total  field  decay  is 
of  the  order  of  50-55  dB.  The  rate  of  decay  of  H^^  computed  by  either 
technique  equals  3.3  dR/ciii,  which  compares  favorably  with  the  3.46  dB/cm 
rate  predicted  by  simple  waveguide  mode  theory  for  this  beyond-cutoff  case. 

Figure  13  graphs  contour  maps  of  the  FD-TD  computed  field  components 
at  the  horizontal  symmetry  plane  of  the  cylinder  1n  an  attempt  to  Indicate 
the  level  of  detail  achieved. 


The  computer  memory  and  running  time  remains  unchanged  If  the  cavity  Is 
loaded  by  aiiiltrary  metal  and  dielectric  contents,  for' a  fixed  number  of 
time  steps. 
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re  10.  FO-TO  MODEL  ffiOWETRY  OF  Figure  11.  TRANSVERSE  CROSS  SECTION 

N-ENDED  ALlfllNUM  CYLINDER  AT  OF  CYLINDER  MODEL  OF  FIGURE  10 

HORIZfflITAL  SYMMETRY  PLANE 


COJffnjTED  HELD  CONTOURS  IN  HORIZWITAL  SYMMETRY  PLANE  OF  CYLINDER  OF  H GORES  10  AND  11 


3.1,2  Loaded  Missile  Guidance  Sectlont  Axial  I ncldence.  Case  [13]. 

This  problem  Involved  the  computation  of  the  fields  within  a  detailed 
model  of  a  metal -coated  missile  guidance  section,  as  shown  1n  cross-section 
in  Figures  14  and  15,  Illuminated  at  axial  Incidence  by  the  300  MHz  plane 


wave,  (E, 


‘Inc  ^1nc 


). 


Two  apertures  were  modeled:  a  circular  one  In  the  nose  (behind  the 
magnesium  fluoride  Infrared  dome);  and  a  circumferential  sleeve  fitting 
23  cm  aft.  The  missile  body  beyond  the  sleeve  fitting  was  assumed  to 
continue  to  Infinity  with  a  constant  cross-section  shape.  The  following 
metal  and  dielectric  Interior  components  were  modeled: 


1)  Head  coll  assembly  (assumed  solid  metal); 

2)  Cooled  detector  unit  -  COU  (assumed  solid  metal); 

3)  Phenolic  ring  around  the  CDU; 

4)  Preamp  can  (metal); 

5)  Wlre^  connecting  the  CDU  to  the  preamp  can; 

6)  Wlre^  connecting  the  preamp  can  to  the  metal  backplane;  and 

7)  Longitudinal  metal  support  rods. 


The  following  electrical 
Ing  the  model : 

Medium 

Aluminum 

Fiberglas 

Phenolic 

Magnesium  fluoride 


parameters  were  assumed 

Relative 

Permittivity, 

1.0 

5.5 

4.5 
5.3 


for  the  media  comprls- 

Conductlvlty 
n  (mhos/ni) 

3.7'10^ 

2.4‘10‘^ 

a.O'lO'^ 

0.0 


The  two  wires  were  really  Idealizations  of  a  more  complex  situation  In 
which  two  mul t1 conductor  wire  bundles  extended  between  the  structures 
mentioned.  By  using  a  simple,  single-wire  model  for  each  bundle,  only 
the  bulk  current  of  each  bundle  (net  metal -to-metal  current)  was  modeled 
using  the  FD-TD  method.  The  wire  model  was  that  of  Section  2. 2. 7.1  of 
this  report. 
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Figure  14.  GEOMETRY  OF  GUIDANCE  SECTION 
MODEL  AT  VERTICAL  SYMMETRY  PLANE. 
SHOWING  COMPONENT  MATERIALS 


phvnoUc  cup 


Figure  15.  GEOMETRY  OF  GUIDANCE 
SECTION  MODEL  AT  HORIZONTAL 
OBSERVATION  PLANE 


Since  only  first-order  correct  radiation  conditions  were  available  at  the 
time,  an  air  conductivity,  a  =  0.025  mho/iti,  was  used  to  attenuate  scattered 
waves.  To  accelerate  convergence  of  the  Interior  fields,  an  air  conductivity, 
0  *  0.025  luho/in,  was  applied  to  the  total  fields  within  the  guidance  section. 

The  model  was  implemented  on  a  24  x  100  x  48-cell  (690,000  unknown  field 
components)  FD-TD  problem  lattice.  A  lattice  resolution  6  »  1/3  cm  =  Aq/300 
was  used,  with  a  time  step  6t  »  6/2c  »  5.55  ps.  1.6  *  10®  words  of  memory  and 
7.0  min  of  Control  Data  Star-lOO  central  processor  time  were  required  to  com¬ 
plete  an  1800  time-step  run  (3.0  cycles  of  the  Incident  wave). 

Figure  16  graphs  contour  maps  of  the  FO-TD  computed  field  components  at 
the  vertical  symmetry  plane  of. the  guidance  section.  An  Important  obser¬ 
vation  Is  that  the  wires  connecting  the  cooled  detector  unit,  preamp  can, 
and  metal  backplane  are  paralleled  by  high-level  magnetic  field  contours 
(Figure  16b).  This  Is  Indicative  of  substantial,  uniform  current  flow  along' 
each  wire.  Such  current  flow  would  generate  locally  a  magnetic  field  looping 
around  the  wire  which,  when  "cut"  by  the  vertical  symmetry  plane,  shows  up 
as  parallel  field  contours  spaced  equally  on  each  side  of  the  wire.  Using 
a  simple  Ampere's  law  argument,  can  be  computed  as  being  equ,il  to 

H‘2iTr  ,  where  H  Is  the  magnitude  of  the  parallel  magnetic  field  contour  and 
r  Is^its  separation  from  the  wire  center.  In  this  manner,  Table  1  lists 
the  predicted  values  of  (current  1n  the  wire  from  the  cooled  detector 
unit  to  the  pre-amp  can)  and  I2  (current  In  the  wire  from  the  pre-amp  can 
to  the  backplane).  For  the  case  of  It  Is  assumed  that  H  -  and 
r  B  i.5  '  6  ■  0.005  m.  For  the  case  of  Ig,  H  is  assumed  that  H  » 

and  r^j  ■  1 .5' 6  =  0.005  m. 
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COMPUTED  FIELD  CONTOURS  IK  VERTICAL  SYMMETRT  PLANE  OF  GUIDANCE  SECTION  OF  FI 


TABLE  r 


PREDICTED  GROUND  WIRE  CURRENTS 


P^n(.(uW/cm2) 

I,()iA) 

IglpA) 

0.1 

2.65-10“^ 

2. 65- 10'^ 

8.3 

26.3 

1.0 

2.65*10'^ 

0.265 

83 

263 

10.0 

2.65’10“^ 

26,5 

030 

2630 

The  inlsflle  guidance  section  model  demonstrates  the  capability  of  the 
pure  TD-TD  method  to  map  fields  penetrating  Into  a  complex  structure  that 
has  multiple  apertures  and  Interior  dielectric  and  metallic  materials.  De¬ 
termination  of  the  accuracy  bound  for  this  method,  as  applied  to  such  com¬ 
plex  structures,  awaits  the  results  of  future  experimental  prog;'.ims  since 
alternative  numerical  approaches  cannot  deal  with  this  level  of  complexity. 

3.2  Hybrid  MOM/FD-TO  Method.  Three-Dimensional  Penetration  Problem-- 
Loaded  Missile  Guidance  Section,  Axial  Incidence  Case  [13],  [34], 
lAsj _ _ _ _ _ 

The  MOM/FD-TD  hybrid  model  for  the  axial -incidence  case  of  the  loaded 
missile  guidance  section  (as  shown  in  Figures  17  and  !4,  15)  was  run  for 
1800  timp  steps  (3.G  cycles  of  the  Incident  wave  at  F  »  300  MHz).  The  ex¬ 
citation  consisted  of  data  for  the  magnitude  and  phase  of  the  electric 
current  distribution  over  the  loci  of  the  short-circuited  nose  aperture  and 
sleeve- fitting  aperture  (as  required  by  Schelkunoff ' s  third  theorem)  com¬ 
puted  using  a  MOM  body-of-revolut1on  code  [37].  The  results  were  compared 
to  pure  FD-TD  results  presented  above  In  an  attempt  to  establish  the  con¬ 
sistency  of  the  MOM/FD-TD  model  for  this  case. 

Figure  Ifi  plots  the  comparison  of  the  hybrid  MOM/FD-TD  results  for  the 
field  contours  In  the  vertical  symmetry  plane  with  the  pure  FD-TD  results 
already  presented  1n  Figure  I6b.  It  Is  seen  that,  for  both  methods  stepped 
to  1800  time  steps,  there  Is  an  excellent  agreement  for  the  0  dB  contours 
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frM  Front  Apertnre  (.»> 


near  the  wire  connecting  the  cooled  detector  unit  to  the  pre-amp  can. 

For  this  contour,  the  maximum  spatial  shift  is  only  about  0.16  cm  (''^0.55) 
in  a  direction  further  out  from  the  wire.  Further,  there  is  observed  to 
be  excellent  agreement  of  the  +10  dB  contour  near  the  wire  connecting  the 
pre-amp  can  to  the  metal  backplane.  For  this  contour,  the  maximum  spatial 
shift  is  only  about  0.1  cm  ('V'0.36)  in  a  direction  further  out  from  the  wire. 
This  implies  that  the  currents  in  these  two  major  wires  are  predicted  to 
be  almost  the  same  by  both  the  pure  FD-TD  method  and  the  hybrid  MOM/FD-TD 
method. 

Figure  19  plots  the  comparison  of  the  hybrid  results  for  the  and 
fields  along  with  a  vertical  cut  tlirough  the  center  of  the  guidance 
section  at  a  point  21  cm  in  back  of  the  nose  aperture  (about  2  cm  in  front 
of  the  sleeve  fitting,  at  the  point  where  the  circumferential  slot  opens 
into  the  interior  of  the  nose  cone).  The  hybrid  run  results  are  after  1500 
time  steps,  while  the  pure  FD-TD  results  are  after  1800  time  steps.  For 
this  case,  a  very  high  level  of  agreement  is  observed  between  the  two  sets 
of  data  at  all  points  of  comparison.  The  worst-case  difference  between  the 
results  is  only  1  dB,  with  most  results  consistent  within  only  fractions  of 
a  decibel 

This  case  study  shows  that  the  hybrid  MOM/FD-TD  approach  yields  results 
for  the  missile  guidance  section  wire  current  and  electromagnetic  fields 
which  are  consistent  with  the  pure  FD-TD  data.  The  study  implies  that  a 
wire  passing  very  close  to  an  aperture,  and  stronnly  coupled  to  that  aper¬ 
ture  (as  for  the  case  of  the  pre-amp  can-to-backpl ane  wire  near  the  sleeve¬ 
fitting  aperture),  can  be  consistently  modeled  using  the  pure  FD-TD  approach 
- - - 

Comparison  of  the  data  sets  for  the  1800  time-step  oaso  for  each  set  results 
in  slightly  lessened  agreement,  such  that  the  worst-case  difference  is  about 
1.7  d3.  This  may  result  from  the  hybrid  program  progressing  to  the  sinusoi¬ 
dal  steady  state  at  a  slightly  faster  rate  than  tne  pure  FD-TD  program,  . 
since  a  sinusoidal  steady-state  equivalent  aperture  current  excitation  is 
employed  from  the  very  beginning  of  the  hybrid  program,  rather  than  an 
aperture  excitation  which  must  build  to  tne  steady  state. 
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and  the  hybrid  MCM/FD-TO  approach.  In  the  latter,  it  must  be  remembered 
that  the  M0M-der1ved  equivalent  aperture  excitation  takes  into  account 
none  of  the  interior  details  of  the  structure.  Finally,  the  study  implies 
that  complex  sub-sections  of  large,  simply-shaped  structures  are  prime 
candidates  for  detailed  modeling  of  the  interior  penetrating  fields  via 
the  hybrid  MOM/FD-TD  method. 

3  Hybrid  FD-TD  Method,  Two-Dimensional  Scattering  Problems 

In  order  to  validate  the  feasibility  of  the  new  formulation  of  the 
FD-TD  lattice  regions  (Section  2.2.2),  the  new  variable-angle  plane  wave 
source  condition  (Section  2.2.2),  the  new  second-order  correct  lattice 
truncation  conditions  (Section  2.2.3),  the  new  sinusoidal  steady-state 
magnitude  and  phase  computation  (Section  2.2.4),  and  the  new  far-field 
scattering  computation  via  the  near-to-far  field  transformation  (Section 
2,2.5),  several  canonical,  two-dimensional,  conducting  and  dielectric 
structures  were  studied  during  the  present  research  program.  The  numerical 
results^  of  the  FD-TD  computed  surface  electric  current  distribution  and 
near  electric  and  magnetic  fields  are  presented  for  the  case  of  a  two- 
dimensional,  square  metal  cylinder  subject  to  TM-polarized  Illumination 
at  normal  and  oblique  incidence.  These  electric  currents  and  near  fields 
are  compared  to  the  MOM-computed  results  [21,  22],  The  scattered-field 
pattern  and  the  corresponding  radar  cross  section  (RCS)  are  derived  from 
the  near-to-far  field  transformation  of  the  FD-TD  data.  These  are  then 
compared  to  the  results  obtained  by  using  the  MOM.  Additional  RCS  results 
for  the  case  of  circular  metal  and  dielectric  cylinders  are  also  presented. 

It  is  shown  that  a  very  high  degree  of  correspondence  is  obtained  using 
this  method. 

■*'*3.3.1  Square  Metal  Cylinder,  Normal  (Broadside)  Incidence, 

TM  Polarization  of  Incident  Wave _ _ 

We  first  consider  the  example  of  the  scattering  of  a  plane  wave  by 
a  square  conducting  cylinder.  The  cylinder  has  the  electrical  size  s  =  2, 
where  s  is  the  width  of  the  side  of  the  cylinder.  The  plane-wave  excitation 


Hhe  numerical  results  were  obtained  using  the  FD-TD  computer  program 
listed  in  Appendix  A.  Standard  Fortran  is  used  for  this  program. 
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is  TM  polarized,  with  field  components  and  !i\  and  propagates  in  the 

z  ^ 

+y  direction  so  that  it  Is  at  normal  Incidence  to  one  side  of  the  cylinder 
(41^  *  0*).  An  84-po1nt  MOM  solution  of  Equation  (18a)  1s  used  as  the  bench¬ 
mark  for  comparison  with  all  FD-TD  results,  with  pulse  current  expansion  and 
point  matching  [21,  22]. 

For  the  FD-TD  analysis,  the  square  cylinder  Is  embedded  in  a  two-dimen¬ 
sional  lattice  as  shown  in  Figure  4.  Each  side  of  the  cylinder  spans  20 
lattice-cell  divisions.  The  connecting  virtual  surface  between  the  FD-TD 
total-field  and  scattered-field  regions  Is  located  at  a  uniform  distance  of 
5  cells  from  the  cylinder  surface.  Figure  20  shows  the  geometry  of  the 
square  cylinder,  the  Incident  plane  wave,  and  the  loci  abed  (cylinder  surface) 
and  ABCD  (off-surface,  near-scattered  field)  along  which  comparative  FD-TD 
and  MOM  results  will  be  graphed. 

Figures  21a  and  21b  graph  the  comparative  results  for  the  FD-TD  and 
MOM  analyses  of  the  magnitude  and  phase  of  the  cylinder  surface  electric 
current  distribution  for  this  case.  Here,  the  FD-TQ  computed  surface  current 
Is  taken  as  n  x  where  n  Is  the  unit  normal  vector  at  the  cylinder  sur¬ 
face,  and  Is  the  FD-TD  value  of  the  magnetic  field  parallel  to  the 
cylinder  surface,  but  at  a  distance  of  0.5  space  cell  from  the  surface.  (The 
displacement  of  the  H  component  from  the  cylinder  surface  Is  a  consequence 
of  the  spatially-interleaved  nature  of  the  E  and  H  components  of  the  FD-TD 
lattice,  Indicated  1n  Figure  2).  In  Figure  21a,  the  magnitude  of  the  FD-TD 
computed  surface  current  agrees  with  the  84-po1nt  MOM  solution  to  better 
than  (0.09  dB)  at  all  comparison  points  more  than  2  cells  from  the 
cylinder  edges  (current  singularities).  In  Figure  21b,  the  phase  of  the 
FD-TD  solution  agrees  with  the  MOM  solution  to  within  +3°  at  all  comparison 
points,  Including  the  shadow  region.  The  uncertainty  bars  shown  1n  this 
figure  Indicate  the  present  level  of  Imprecision  1n  using  the  FD-TH  method 

.y 

to  locate  the  constant-phase  points  of  the  computed  time-domain  wave¬ 
form  (equivalent  to  +1  time  step).  This  Imprecision  can  be  reduced  in  future 
FD-TD  programs  by  Incorporating  a  simple  Interpolation  algorithm  to  achieve 
fractional  time-step  resolution  of  points  of  constant  phase. 

Figures  22a,  22b,  23a  and  23b  show  the  comparison  of  the  magnitude  and 
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Fig.  23b.  COMPARISON  OF  MOM  AND  FD-TD  RESULTS 
FOR  PHASE  OF  NEAR  MAGNETIC  FIELDS 
TANGENTIAL  TO  CONTOUR  So 


phase  of  the  near  scattered  electric  and  magnetic  fields  computed  by  the 
FD-TD  method  and  MOM.  The  electric  field  is  tangential  to  a  virtual 
surface,  located  at  a  uniform  distance  of  h  ■  7  space  cells  from  the 
cylinder  surface;  and  the  magnetic  field  Is  tangential  to  a  virtual  surface, 
S'^,  located  at  6.5  cells  from  the  cylinder  surface.  Both  virtual  surfaces 
are  embedded  In  the  scattered-field  region  of  the  FD-TD  lattice.  The  level 
of  correspondence  between  the  FD-TD  and  MOM  results  Is  +2.5%  (+0.2  dB)  and 
+3° . 

In  order  to  obtain  the  far-fleld  pattern  and  the  scattering  cross- 
section,  the  near-field  to  far-fleld  transformation  discussed  In  Section 
2.2.5  Is  followed  for  two-dimensional  structures.  The  near  scattered 
fields  shown  In  Figures  22a,  22b,  23a  and  23b  are  converted 'into  the  corre¬ 
sponding  equivalent  electric  and  magnetic  current  distributions  along 
according  to  Equations  (14a)  and  (14bl.  The  far  field  Is  then  computed 
by  applying  Equations  (17a)  -  (17e). 

Figure  24  plots  the  normalized  radar  cross  section  (RCS)  versus 
scattering  angle,  41,  for  the  square,  conducting  cylinder,  Two  cases  aro 
plotted.  First,  the  solid  curve  represents  the  RCS  computed  directly  from 
the  Induced  surface  electric  currents  of  the  evlinder.  These  surface 
currents  were  derived  via  a  MOM  solution  of  Equ^tljn  (18),  and  the  far  field 
derived  by  performing  the  Integration  of  Equatloii  (19).  Second,  the  dots 
represent  RCS  values  computed  by  converting  the  KD- 'D  derived  near  electric 
and  magnetic  fields  of  Figures  22a,  22b,  23a,  23b  to  the  far  field  by 
applying  Equations  (14a),  (14b),  end  (17a)  -  (17e).  A  very  high  degree  of 
correspondence  1s  Indicated  between  the  two  sets  of  results.  This  corre¬ 
spondence  Indicates  the  validity  of  the  hybrid  FD-rn  method  for  computing 
the  far  flex'd  and  tadar  cross  section,  as  well  as  the  accuracy  of  the  new 
FD-TD  algorithm  provisions. 
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**3.3.2  Square  Metal  Cylinder,  Oblique  (45“)  Incidence,  TM 
Polarization  of  Incident  Wave _ _ 

We  next  consider  the  oblique  illumination  of  the  square  metal  cylinder 
of  the  previous  exampl®.  Figure  25  shows  the  geometry  of  the  square  cylinder 
the  incident  plane  wave,  and  the  loci  (cylinder  surface)  and  WC  (off- 
surface,  near-scattered  field)  along  which  comparative  FD-TD  and  MOM  results 
win  be  graphed.  The  only  change  made  from  the  previous  scattering  example 
is  that  the  incident  wave  is  assumed  to  propagate  at  an  angle  of  45°  rela¬ 
tive  to  one  side  of  the  cylinder  (41^  '  45°).  This  change  is  made  simply  by 
altering  one  data  statement  of  the  Fortran  computer  program. 

Figures  26a  and  26b  graph  the  comparative  results  for  the  FD-TD  and 
MOM  analyses  of  the  magnitude  and  phase  of  the  cylinder  surface  electric 
current  distribution  for  this  case.  Again,  the  FD-TD  computed  surface 
current  is  taken  as  n  x  where  n  is  the  unit  normal  vector  at  the 
cylinder  surface,  and  is  the  FD-TD  value  of  the  magnetic  field  parallel 
to  the  cylinder  surface,  but  at  a  distance  of  0.5  space  cell  from  the  surface 
In  Figure  26a,  the  magnitude  of  the  FD-TD  computed  surface  current  agrees 
again  with  the  84-po1nt  MOM  solution  to  better  than  +I5i  (+0.09  dB)  at  all 
comparison  points  more  than  2  cel  Is  from  the  cylinder  edges  (current  singu¬ 
larities).  In  Figure  26b,  the  phase  of  the  FD-TD  solution  again  agrees  with 
the  MOM  solution  to  within  +3°  at  virtually  each  comparison  point,  including 
the  shadow  region. 

Figure  27a,  27b,  28a,  and  28b  show  the  comparison  of  the  magnitude 
and  phase  of  the  near-scattered  electric  and  magnetic  fields  co*-  uted  by  the 
FD-TD  method  and  MOM.  Again,  the  electric  field  is  tangential  to  a  virtual 
surface,  S^,  located  at  a  uniform  distance  of  h  7  space  cells  From  the 
cylinder  surface;  and  the  magnetic  field  is  tangential  to  a  virtual  surface, 
S-*g,  located  at  6.5  cells  from  the  cylinder  surface.  The  level  of  corre- 
pondencB  between  the  FD-TD  and  MOM  results  is  again  ±Z,5%  (+0.2  dB)  and  +3°. 
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GEOMETRY  FOR  TWO-DIMENSIONAL  SCATTERING 
EXAMPLE  :  SQUARE  CONDUCTING  CYLINDER 
ILLUMINATED  BY  A  PLANE  WAVE  AT  OBLIQUE 
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Fig.  26a.  COMPARISON  OF  MOM  AND  FD-TD  RESULTS 
FOR  MAGNITUDE  OF  ELECTRIC  CURRENTS  ON 
SURFACE  OF  CYLINDER,  OBLIQUE  INCIDENCE 
CASE  ((^  =  45») 
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Fig  26b.  COMPARISON  OF  MOM  AND  FD-TO  RESULTS 
FOR  PHASE  OF  ELECTRIC  CURRENTS  ON 
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Fig.  27a.  COMPARISON  OF  MOM  AND  FD-TD  RESULTS 
FOR  NEAR  ELECTRIC  FIELD  TANGENTIAL  TO 
CONTOUR  Sq  for  OBLIQUE  INCIDENCE 
CASE  {(p  ‘  45") 

72 


« 


_  mom  (80-Point  Solution  ) 

^  ^  ^  FO-TD  (6-Cycle  Solution) 


Position  On  Contour  Sq 

)mparison  of  mom  and  fd -to  results 

)R  PHASE  OF  NEAR  ic 

\NGENTIAL  TO  CONTOUR  Sq  FOR  OBLIQUE 

nnFwrF  CASE  ( cfc  =  45®) 


-  MOM  (80- Point  Solution) 

+  PT-TD(6-Cycle  Solution) 


Position  On  Contour  Sq 

Fig.  28b.  COMPARISON  OF  MOM  AND  FD-TD  RESULTS 
FOR  PHASE  OF  NEAR  MAGNETIC  FIELD 
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The  high  level  of  agr'eement  between  the  FD-TD  and  MOM  results  for 
both  the  =  and  4'^'’ 45°  angles  of  incidence  verifies  that  the  key,  new 
FD-TD  features  (variable-incidence  plane  wave,  second-order-correct  lattice 
truncations,  sinusoidal  steady-state  magnitude  and  phase  computations)  work, 
and  work  very  well.  In  fact,  the  observed  level  of  accuracy  for  these 
latest  results  is  as  much  as  5  to  10  times  as  good  as  the  previously-observed 
[12],  [13]  FD-TD  accuracy  level  of  +10%  (+1  dB).  It  is  expected  that  the 
new,  enhanced  accuracy  levels  will  be  observed  when  the  present  two-dimen¬ 
sional  FD-TD  computer  programs  are  modified  to  the  three-dimensional  case. 

**3.3.3  Verification  of  Near-to-Far  Field  Transformation  for 
Circular  Metal  and  Dielectric  Cylinders  _ _ _ 

To  provide  further  verification  of  the  application  of  the  near-to-far 
field  transformation  of  Section  2.2.5,  pure  MOM  studies  have  been  performed 
during  this  research  program  on  a  metal  and  a  dielectric  circular  cylinder. 
These  studies  concern  the  computation  of  scattering  from  such  cylinders  for 
the  two-dimensional  case  of  plane  wave  Illumination  (TM  polarization). 

Figure  29  shows  the  normalized  radar  cross  section  (RCS)  computed  for 
a  metal  cylinder  of  size  kga  »  1 .  The  solid  line  shows  the  RCS  computed 
directly  from  the  surface  electric  current,  J21  on  the  cylinder.  The  dots 
show  the  RCS  computed  using  the  near-to-far  field  transformation  executed 
along  a  square  virtual  surface  located  at  an  electrical  distance  k^h  =  0.38 
from  the  surface  of  the  cylinder.  The  agreement  of  the  two  solutions  for 
the  RCS  is  extremely  high  (to  better  than  4  decimal  places  at  each  point  of 
comparison) . 

Figure  30  shows  the  normalized  RCS  computed  for  a  dielectric  cylinder 
of  size  k^a  =  .63.  The  solid  line  shows  the  RCS  computed  directly  from 
the  equivalent  electric  and  magnetic  currents  on  the  surface  of  the  cylinder, 
which  had  been  previously  derived  using  a  MOM  solution  of  a  surface  integral 
equation.  The  dots  show  the  RCS  computed  using  the  near-to-far  field  trans¬ 
formation  executed  along  a  square  virtual  surface  located  at  an  electrical 
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Figure  30.  EQUIVALENCE  OF  COMPUTED  RESULTS  FOR  MDAR  CROSS  SECTION  OBTAINED  USING  THE  HOH  SURFACE 
CURRENTS  APPROACH  .AND  THE  HYBRID  MOM  HEAR-TO-FAR  FIELD  TRANSFORMATION  TECHNIQUE 


distance  k^h  ®  0.24  from  the  surface  of  the  cylinder.  Again,  the  agreement 
of  the  two  solutions  for  the  RCS  1s  extremely  high  (to  better  than  4  decimal 
places  at  each  point  of  comparison). 

It  Is  therefore  concluded  that  the  application  of  the  near-to-far 
field  transformation  along  a  rectangular  locus  surrounding  a  scatterer  Is, 
for  all  practical  purposes,  an  exact  procedure  for  deriving  the  far  fields. 
The  validity  of  this  near-to-far  field  transformation  Is  the  heart  of  the 
proposed  hybrid  FD-TD  method  for  computing  the  far  scattered  fields  and  RCS 
from  arbitrary,  extremely-complex  structures  in  three  dimensions.  It  should 
be  noted  that  this  procedure  should  also  be  equally  valid  and  useful  In 
computing  the  far-fleld  radiation  patterns  of, extremely-complex  source 
regions  and  antennas. 
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4.0  GUMMARY  AND  CONCLUSIONS 

Electromagnetic  penetration  and  scattering  problems  are  difficult  to 
treat  with  many  analytical  or  numerical  methods  because  of  the  Inability  of 
these  methods  to  simply  deal  with  the  effects  of  structure  materials,  aper¬ 
tures,  curvatures,  corners,  and  internal  contents.  In  two  earlier  RADC 
contracts,  IITRI  Investigated  the  application  of  a  new  approach  for  the 
direct  modeling  of  very  complex  electromagnetic  interaction  problems:  the 
finite-difference,  time-domain  (FD-TD)  solution  of  Maxwell's  equations. 

The  FD-TD  method  has  key  advantages  relative  to  available  model ing  approaches. 
These  advantages  permit  it  to  accurately  treat  complex  problems  that  are 
beyond  the  scope  of  solution  by  any  other  method.  The  ultimate  aim  of 
research  in  this  area  is  to  develop  an  accurate,  easily-used,  general  com¬ 
puter  program  solving  for  either  electromagnetic  field  penetration,  scatter¬ 
ing,  or  radiation  for  arbitrary  metal/dielectric  structures  spanning  up  to 
10  nr  more  wavelengths  in  three  dimensions  with  a  spatial  resolution  better 
than  0.1  wavelength. 

In  order  to  more  fully  determine  the  usefulness  of  the  FD-TD  method, 

RADC  thought  that  it  is  desirable  to  distribute  this  technique  to  as  wide 
a  range  of  users  as  possible  so  that  it  can  be  tested  by  actual  implementa¬ 
tion.  The  overall  objectives  of  algorithm  development  in  this  case  are  to 
allow  RADC  to  write  a  user-oriented  computer  program  for  the  FD-TD  technique. 

The  goals  of  the  present  IITRI  research  effort  for  RADC  include  the 
development  of  specific  algorithms  of  high  importance  to  help  provide  a 
flexible,  s1mple-to-use,  and  highly  accurate  user-oriented  FD-TD  computer 
program.  To  moet  these  goals,  IITRI  has  tested  five  key  improvements  in 
the  FD-TD  algorithm  during  this  effort,  including  the  following: 

1  ■  Total -field/ scattered-field  lattice  di vi sion . 

This  permits  a  very  high  computational  dynamic  range  to  ac¬ 
curately  model  fields  within  shadow  zones  or  cavities.  This 
further  permits  programming  of  variable  angle  of  incidence 
and  the  second-order  correct  radiation  condition,  summarized 
below. 
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Variable  ajacj  1  e  of  incidence. 

For  two-dimensional  problems,  this  permits  a  single  data  card  or 
Fortran  statement  to  specify  1n  a  very  accurate  manner  the  angle 
of  incidence  of  a  plane  wave  lllutnlnating  a  structure.  For  three- 
dimensional  problems,  both  the  angle  of  Incidence  and  polarization 
could  be  specified.  There  Is  no  requirement  to  rotate  the  geometry 
of  the  Interacting  structure  In  the  FD-TD  lattice. 

3 •  Second-order  accurate  radiation  condition . 

This  reduces  the  uncertainty  of  the  final  computed  results  by  as 
much  as  ten-to-one.  FD-TD  computations  using  this  radiation  con¬ 
dition  now  have  estimated  field-magnitude  uncertainties  of  better 
than  +2.5',ll  (+0.2  dB)  versus  previous  uncertainties  of  _+10'X-_+15‘y. 

(+.1  dB). 

4.  Magnitude  and  phase  computation  condition  for  the  sinusoidal  steady 

state,. _ _ _ _ _ _ _ _ _ 

This  permits  accurate  determination  of  the  magnitude  and  phase  of 
FD-TD  computed  fields  at  any  desired  points  for  later  use  In  com¬ 
putations  involving  scattering,  radiation  or  coupling  to  wires. 

This  approach  avoids  any  ambiguity  due  to  either  a  possible  DC 
offset  of  the  fields  or  the  repetitive  nature  of  the  sinusoidal 
waveform. 

5 .  Nea  r-  to  -far  _f  i,e  1  d  ,t  r  a  n s  formation. 

This  permits  the  far  scattered  fields  and  radar  cross  section  of 
arbitrary  structures  modeled  by  the  FD-TD  method  to  be  easily  and 
accurately  determined.  Observed  accuracy  of  the  radar  cross  section 
using  this  feature  Is  In  the  order  of  +1X  (+0.09  dB). 

In  addition  to  the  five  algorithm  Improvements  tested  byllTRI,  this  report 
also  summarized  a  FD-TD  feature  which  has  recently  appeared  In  the  literature 
that  permits  computation  of  the  coupling  of  currents  to  thin  wires  and  struts. 
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The  conclusions  of  this  report  .'.re  ds  follows: 

1.  The  accuracy  of  the  pure  FD-TD  method  for  electromagnetic  Interaction 
problems  can  reach  the  high  levels  previously  attained  only  by  method- 
of-moments  (MOM)  approaches  when  the  second-order  accurate  radiation 
condition  is  used  in  the  FD-TD  algorithm.  The  FD-TD  method  retains  its 
key  advantages  over  MOM  in  terms  of  the  much  larger  electrical  size 
and  greater  complexity  of  the  structures  that  can  be  modeled. 

2.  The  specification  of  variable  angle  of  Incidence  and  polarization  of 
an  illuminating  wave  can  he  achieved  with  the  FD-TD  method  using  only 
a  single  data  card  or  Fortran  statement.^ 

3.  The  total-field/scattered-field  regional  division  of  the  FD-TD  lattice 
can  be  successfully  Implemented  and  offers  the  significant  advantage 
of  a  high  computational  dynamic  range.  In  additioni  this  lattice  divi¬ 
sion  provides  a  framework  for  programming  variable  wave  incidence  and 
polarization,  improved  radiation  conditions,  and  the  near-field  to 
far-field  transformation  for  scattering  problems. 

4.  The  near-to-far  field  transformation  along  a  'rcctangul ar  virtual  sur¬ 
face  surrounding  a  scatterer  makes  it  possible  to  use  the  FD-TD  method 
to  compute  the  far  scattered  fields  and  radar  cross  section  of  complex, 
arbitrary  structures  with  great  precision. 

It  is  the  opinion  of  the  authors  of  this  repo>"t  that  the  FD-TD  method  deserves 
additional  investigation  to  probe  just  what  are  the  limits  of  application  of 
this  ext*'3niely  prrr.  '  ng  approach  to  accurately  model  electromagnetic  penetra¬ 
tion,  scattering,  and  radiation  problems. 


This  incident  wave  specification  is  now  as  simple  for  the  FD-TD  method  as  it 
has  been  with  MOM.  Howevor,  the  FD-TD  approach  requires  re-running  the  entire 
problem  for  each  new  incident  wave  angle.  With  MOM,  only  a  single  inversion 
of  the  system  matrix  is  required.  Subsequently,  arbitrary  wave  excitation  is 
treated  as  a  simple  matrix  multiplication  of  the  equivalent  excitation 
vector.  MOM  therefore  permits  a  conceptually  simple'"  treatment  of  the  vari¬ 
able  wave  incidence  problem. 
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APPENDIX  A 


Listing  of  Fortran  Computer  Program 


This  appendix  contains  the  listing  of  the  Fortran  computer  program  used 
to  obtain  the  FD-TD  results  for  the  square-cylinder  scattering  problems  of 
Section  3.3  of  this  report.  This  computer  program  contains  all  of  the  new 
FO-TD  algorithms  that  have  been  examined  during  the  present  research  effort, 
except  the  thin-wire  coupling  algorithm  of  Section  2. 2. 7. 2.  This  program 
provides  a  general  treatment  of  two-dimensional  penetration  and  scattering 
problems  for  the  transverse  magnetic  (TM)  polarization  case,  and  utilizes  a 
50  x  49-cell  Yee  grid.  The  following  is  a  succinct  description  of  the  user 
Inputs  to  this  program. 


Program  Line 

Symbol  or  Data  Input 

Function 

004 

FREQ 

Illumination  frequency  (Hz) 

005 

DX 

Lattice  cell  size  (meters) 

006 

NCYCS 

Number  of  complete  cycles  of 
the  sinusoidal  incident  wave 
that  the  program  is  time-stepped 

007 

MEDIA 

Equal  to  2  +  number  of  distinct 
dielectric  or  conducting  mater¬ 
ials  modeled  in  the  FD-TD  grid 

008 

DATA  EPS 

Relative  permittivity  of  mater¬ 
ials  modeled.  User  specification 
begins  wUh  data  item  #4  in 
list.  External  air  is  item  #1. 

009 

DATA  SIG 

Electric  conductivity  (mhos/m) 
of  materials  modeled.  User  spec¬ 
ification  begins  with  data  item 
#4  in  list.  Air  is  item  #1 . 

OlO 

PHI 

Angle  of  incidence  (degrees) 
relative  to  y  axis  of  grid.  Lim¬ 
ited  to  range  0"  -  90“  in  pres¬ 
ent  program. 

on 

ISA 

FO-TD  "1"  coordinate  of  left 

field  matching  plane  (between 
total-field  and  scattered-field 
regions) . 
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Program  line 


Symbol  or  Data  Input 


Function 


103  -  10/ 


MEZ(I,J) 


FD-TD  ‘M"  coordinate  of  right 
field-matching  plane  (between 
total -field  and  scattered-field 
regions) . 

FD-TO  "j"  coordinate  of  front 
field-matching  plane  (between 
total-field  and  scattered-field 
regions) . 

FD-TD  "j"  coordinate  of  back 
field-matching  plane  (between 
total-field  and  scattered-field 
regions) , 

Assigns  a  material -medium  type 
number  to  the  field  compo¬ 
nents  of  the  grid  that  comprise 
the  structure  to  be  modeled. 

The  type  number  corresponds  to 
the  EPS  and  SI6  specifications 
of  program  l^nes  DOS  and  009, 
Arbitrary  specification  is  pos¬ 
sible  simply  by  using  any  number 
of  Fortran  statements  here  to 
make  the  assignment.  In  the  ex¬ 
ample  shown,  a  single  DO  loop 
suffices  to  specify  that  the 
object  is  square  and  has  a  sur¬ 
face  composed  of  material -medium 
type  H  (e„  =  1.0,  a  =  3.72  - 
7  " 

10  mhos/meter).  The  side  of  the 
square  object  is  specified  to  be 
20  grid  cells  across  with  cor¬ 
ners  at  (16,16),  (36,16), 

(16,36) ,  and  (36,36) . 
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